DOFVector.hh 39.9 KB
Newer Older
1
#include <list>
Thomas Witkowski's avatar
Thomas Witkowski committed
2
#include <algorithm>
3
#include <math.h>
4
5
6
7
8
9
10
11
12
13
14
15
16
17

#include "FixVec.h"
#include "Boundary.h"
#include "DOFAdmin.h"
#include "ElInfo.h"
#include "Error.h"
#include "FiniteElemSpace.h"
#include "Global.h"
#include "Mesh.h"
#include "Quadrature.h"
#include "AbstractFunction.h"
#include "BoundaryManager.h"
#include "ElementVector.h"
#include "Assembler.h"
18
#include "OpenMP.h"
19
#include "Operator.h"
20
#include "Parameters.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
21
#include "Traverse.h"
22
23
24
25

namespace AMDiS {

  template<typename T>
26
  DOFVectorBase<T>::DOFVectorBase(const FiniteElemSpace *f, std::string n)
Thomas Witkowski's avatar
Thomas Witkowski committed
27
28
29
    : feSpace(f),
      name(n),
      boundaryManager(NULL)
30
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
31
    nBasFcts = feSpace->getBasisFcts()->getNumber();
Thomas Witkowski's avatar
Thomas Witkowski committed
32
    dim = feSpace->getMesh()->getDim();
33

34
    localIndices.resize(omp_get_overall_max_threads());
Thomas Witkowski's avatar
Thomas Witkowski committed
35
    localVectors.resize(omp_get_overall_max_threads());
36
    grdPhis.resize(omp_get_overall_max_threads());
Thomas Witkowski's avatar
Thomas Witkowski committed
37
    grdTmp.resize(omp_get_overall_max_threads());
38
    D2Phis.resize(omp_get_overall_max_threads());
39

40
    for (int i = 0; i < omp_get_overall_max_threads(); i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
41
42
      localIndices[i] = GET_MEMORY(DegreeOfFreedom, this->nBasFcts);
      localVectors[i] = GET_MEMORY(T, this->nBasFcts);
43
      grdPhis[i] = NEW DimVec<double>(dim, DEFAULT_VALUE, 0.0);
Thomas Witkowski's avatar
Thomas Witkowski committed
44
      grdTmp[i] = NEW DimVec<double>(dim, DEFAULT_VALUE, 0.0);
45
      D2Phis[i] = NEW DimMat<double>(dim, NO_INIT);
46
    }
47
48

    elementVector = NEW ElementVector(this->nBasFcts);
49
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
50
  
51
  template<typename T>
Thomas Witkowski's avatar
Thomas Witkowski committed
52
53
54
55
  DOFVectorBase<T>::~DOFVectorBase()
  {
    for (int i = 0; i < static_cast<int>(localIndices.size()); i++) {
      FREE_MEMORY(localIndices[i], DegreeOfFreedom, this->nBasFcts);
Thomas Witkowski's avatar
Thomas Witkowski committed
56
      FREE_MEMORY(localVectors[i], T, this->nBasFcts);
57
      DELETE grdPhis[i];
Thomas Witkowski's avatar
Thomas Witkowski committed
58
      DELETE grdTmp[i];
59
      DELETE D2Phis[i];
Thomas Witkowski's avatar
Thomas Witkowski committed
60
61
62
63
    }
  }
  
  template<typename T>
64
  DOFVector<T>::DOFVector(const FiniteElemSpace* f, std::string n)
65
    : DOFVectorBase<T>(f, n),
66
      coarsenOperation(NO_OPERATION)
67
68
  {
    init(f, n);
69
70
71
72

#ifdef HAVE_PARALLEL_AMDIS
    applicationOrdering = NULL;
#endif
73
74
75
  } 

  template<typename T>
76
  void DOFVector<T>::init(const FiniteElemSpace* f, std::string n)
77
78
79
  {
    name = n;
    feSpace = f;
80
    if (feSpace && feSpace->getAdmin()) {
81
      (feSpace->getAdmin())->addDOFIndexed(this);
Thomas Witkowski's avatar
Thomas Witkowski committed
82
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
83
    this->boundaryManager = NEW BoundaryManager(f);
84
85
86
87
88
  }

  template<typename T>
  DOFVector<T>::~DOFVector()
  {
89
    if (feSpace && feSpace->getAdmin()) {
90
91
92
93
94
95
      (feSpace->getAdmin())->removeDOFIndexed(this);
    }

    if (this->boundaryManager) {
      DELETE this->boundaryManager;
    }
96
97

    vec.clear();
98
99
100
101
102
103
104
105
106
107
108
  }

  template<typename T>
  DOFVector<T> * DOFVector<T>::traverseVector = NULL;

  template<typename T>
  FastQuadrature *DOFVector<T>::quad_fast = NULL;

  template<typename T>
  double DOFVector<T>::norm = 0.0;

Thomas Witkowski's avatar
Thomas Witkowski committed
109
110
111
112
113
114
115
116
117
118
119
120
121
122
  template<typename T>
  void DOFVectorBase<T>::addElementVector(T factor, 
					  const ElementVector &elVec, 
					  const BoundaryType *bound,
					  bool add)
  {
    FUNCNAME("DOFVector::addElementVector()");

    int n_row = elVec.getSize();

    for (DegreeOfFreedom i = 0; i < n_row; i++) {
      BoundaryCondition *condition = 
	bound ? this->getBoundaryManager()->getBoundaryCondition(bound[i]) : NULL;

Thomas Witkowski's avatar
Thomas Witkowski committed
123
      if (!(condition && condition->isDirichlet())) {
Thomas Witkowski's avatar
Thomas Witkowski committed
124
125
126
127
128
129
130
	DegreeOfFreedom irow = elVec.dofIndices[i];
	(*this)[irow] = (add ? (*this)[irow] : 0.0);
	(*this)[irow] += factor * elVec[i];
      }
    }
  }

131
132
133
134
135
  template<typename T>
  double DOFVector<T>::nrm2() const
  {
    FUNCNAME("DOFVector<T>::nrm2()");

136
137
138
    TEST_EXIT_DBG(feSpace)("feSpace is NULL\n");
    TEST_EXIT_DBG(feSpace->getAdmin())("admin is NULL\n");
    TEST_EXIT_DBG(static_cast<int>(vec.size()) >= feSpace->getAdmin()->getUsedSize())
139
      ("size = %d too small: admin->sizeUsed = %d\n", vec.size(),
140
141
142
       feSpace->getAdmin()->getUsedSize());
    
    double nrm = 0.0;
143
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
144
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
145
146
147
148
149
150
151
152
153
154
      nrm += (*vecIterator) * (*vecIterator);

    return(sqrt(nrm));
  }

  template<typename T>
  double DOFVector<T>::squareNrm2() const
  {
    FUNCNAME("DOFVector<T>::nrm2()");

155
156
157
    TEST_EXIT_DBG(feSpace)("feSpace is NULL\n");
    TEST_EXIT_DBG(feSpace->getAdmin())("admin is NULL\n");
    TEST_EXIT_DBG(static_cast<int>(vec.size()) >= feSpace->getAdmin()->getUsedSize())
158
      ("size = %d too small: admin->sizeUsed = %d\n", vec.size(),
159
160
161
       feSpace->getAdmin()->getUsedSize());
    
    double nrm = 0.0;
162
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
163
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
164
165
166
167
168
169
170
171
      nrm += (*vecIterator) * (*vecIterator);

    return nrm;
  }

  template<typename T>
  T DOFVector<T>::asum() const
  {
172
    FUNCNAME("DOFVector<T>::asum()");
173

174
175
176
177
178
    TEST_EXIT_DBG(feSpace)("feSpace is NULL\n");
    TEST_EXIT_DBG(feSpace->getAdmin())("admin is NULL\n");
    TEST_EXIT_DBG(static_cast<int>(vec.size()) >= feSpace->getAdmin()->getUsedSize())
      ("size = %d too small: admin->sizeUsed = %d\n", vec.size(),
       feSpace->getAdmin()->getUsedSize());
179

180
    double nrm = 0.0;
181
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
182
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
183
184
185
186
187
188
189
190
      nrm += abs(*vecIterator);

    return(nrm);
  }

  template<typename T>
  T DOFVector<T>::sum() const
  {
191
    FUNCNAME("DOFVector<T>::sum()");
192

193
194
195
196
197
    TEST_EXIT_DBG(feSpace)("feSpace is NULL\n");
    TEST_EXIT_DBG(feSpace->getAdmin())("admin is NULL\n");
    TEST_EXIT_DBG(static_cast<int>(vec.size()) >= feSpace->getAdmin()->getUsedSize())
      ("size = %d too small: admin->sizeUsed = %d\n", vec.size(),
       feSpace->getAdmin()->getUsedSize());
198

199
    double nrm = 0.0;    
200
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
201
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator)
202
203
204
205
206
207
208
209
      nrm += *vecIterator;

    return(nrm);
  }

  template<typename T>
  void DOFVector<T>::set(T alpha)
  {
210
    FUNCNAME("DOFVector<T>::set()");
211

212
213
214
    TEST_EXIT_DBG(feSpace)("feSpace is NULL\n");
    TEST_EXIT_DBG(feSpace->getAdmin())("admin is NULL\n");
    TEST_EXIT_DBG(static_cast<int>(vec.size()) >= feSpace->getAdmin()->getUsedSize())
215
      ("size = %d too small: admin->sizeUsed = %d\n", vec.size(),
216
217
       feSpace->getAdmin()->getUsedSize());
    
218
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(this), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
219
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator) {
220
      *vecIterator = alpha ; 
Thomas Witkowski's avatar
Thomas Witkowski committed
221
    }
222
223
224
225
226
227
  }


  template<typename T>
  void DOFVector<T>::copy(const DOFVector<T>& x)
  {
228
229
230
231
232
    FUNCNAME("DOFVector<T>::copy()");
    
    TEST_EXIT_DBG(feSpace && x.feSpace)
      ("feSpace is NULL: %8X, %8X\n", feSpace, x.feSpace);
    TEST_EXIT_DBG(feSpace->getAdmin() && (feSpace->getAdmin() == x.feSpace->getAdmin()))
233
234
      ("no admin or different admins: %8X, %8X\n",
       feSpace->getAdmin(), x.feSpace->getAdmin());
235
    TEST_EXIT_DBG(static_cast<int>(vec.size()) >= feSpace->getAdmin()->getUsedSize())
236
      ("size = %d too small: admin->sizeUsed = %d\n", vec.size(), 
237
238
       feSpace->getAdmin()->getUsedSize());
    TEST_EXIT_DBG(static_cast<int>(x.vec.size()) >= feSpace->getAdmin()->getUsedSize())
239
      ("x.size = %d too small: admin->sizeUsed = %d\n", x.vec.size(), 
240
241
       feSpace->getAdmin()->getUsedSize());
    
242
243
    Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(this), USED_DOFS);
    Iterator xIterator(dynamic_cast<DOFVector<T>*>(const_cast<DOFVector<T>*>(&x)), USED_DOFS);
244
245
246
247
248
    for (vecIterator.reset(), xIterator.reset();
	 !vecIterator.end();
	 ++vecIterator, ++xIterator) {
      
      *vecIterator = *xIterator; 
Thomas Witkowski's avatar
Thomas Witkowski committed
249
    }
250
251
252
253
254
255
  }


  template<typename T>
  T DOFVector<T>::min() const
  {
256
257
258
259
260
    FUNCNAME("DOFVector<T>::min()");
    
    TEST_EXIT_DBG(feSpace && feSpace->getAdmin())
      ("pointer is NULL: %8X, %8X\n", this, feSpace->getAdmin());
    TEST_EXIT_DBG((static_cast<int>(vec.size())) >= feSpace->getAdmin()->getUsedSize())
261
      ("size = %d too small: admin->sizeUsed = %d\n", vec.size(), 
262
263
264
       feSpace->getAdmin()->getUsedSize());

    T m;    
265
    Iterator vecIterator(const_cast<DOFIndexed<T>*>(dynamic_cast<const DOFIndexed<T>*>(this)), USED_DOFS);
266
    for (vecIterator.reset(), m = *vecIterator; !vecIterator.end(); ++vecIterator) {
267
      m = std::min(m, *vecIterator);
268
269
    }

270
271
272
273
274
275
    return m;
  }

  template<typename T>
  T DOFVector<T>::max() const 
  {
276
277
278
279
280
    FUNCNAME("DOFVector<T>::max()");
    
    TEST_EXIT_DBG(feSpace && feSpace->getAdmin())
      ("pointer is NULL: %8X, %8X\n", this, feSpace->getAdmin());
    TEST_EXIT_DBG((static_cast<int>(vec.size())) >= feSpace->getAdmin()->getUsedSize())
281
      ("size = %d too small: admin->sizeUsed = %d\n", vec.size(), 
282
       feSpace->getAdmin()->getUsedSize());
283

284
    T m;    
285
    Iterator vecIterator(const_cast<DOFIndexed<T>*>(dynamic_cast<const DOFIndexed<T>*>(this)), USED_DOFS);
286
    for (vecIterator.reset(), m = *vecIterator; !vecIterator.end(); ++vecIterator) {
287
      m = std::max(m, *vecIterator);
288
289
    }

290
291
292
    return m;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
293
294
295
296
297
298
  template<typename T>
  T DOFVector<T>::absMax() const
  {
    return std::max(abs(max()), abs(min()));
  }

299
300
301
302
303
  template<typename T>
  void gemv(MatrixTranspose transpose, T alpha,
	    const DOFMatrix& a, const DOFVector<T>& x,
	    T beta, DOFVector<T>& y)
  {
304
305
306
307
    FUNCNAME("DOFVector<T>::gemv()");

    int j, jcol, ysize;
    T sum, ax;
308
    const DOFMatrix::MatrixRow *row;
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
    
    TEST_EXIT_DBG(a.getRowFESpace() && 
		  a.getColFESpace() && 
		  x.getFESpace() && 
		  y.getFESpace())
      ("feSpace is NULL: %8X, %8X, %8X, %8X\n", 
       a.getRowFESpace(), a.getColFESpace(), x.getFESpace(), y.getFESpace());
    TEST_EXIT_DBG(a.getRowFESpace()->getAdmin() && 
		  a.getColFESpace()->getAdmin() && 
		  (((transpose == NoTranspose) && 
		    (a.getColFESpace()->getAdmin() == x.getFESpace()->getAdmin()) && 
		    (a.getRowFESpace()->getAdmin() == y.getFESpace()->getAdmin()))||
		   ((transpose == Transpose) && 
		    (a.getRowFESpace()->getAdmin() == x.getFESpace()->getAdmin()) && 
		    (a.getColFESpace()->getAdmin() == y.getFESpace()->getAdmin()))))
324
325
326
      ("no admin or different admins: %8X, %8X, %8X, %8X\n",
       a.getRowFESpace()->getAdmin(), a.getColFESpace()->getAdmin(), 
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
327
    
328
    if (transpose == NoTranspose) {
329
      TEST_EXIT_DBG(static_cast<int>(x.getSize()) >= a.getColFESpace()->getAdmin()->getUsedSize())
330
	("x.size = %d too small: admin->sizeUsed = %d\n",
331
332
	 x.getSize(), a.getColFESpace()->getAdmin()->getUsedSize());
      TEST_EXIT_DBG(static_cast<int>(y.getSize()) >= a.getRowFESpace()->getAdmin()->getUsedSize())
333
	("y.size = %d too small: admin->sizeUsed = %d\n",
334
335
	 y.getSize(), a.getRowFESpace()->getAdmin()->getUsedSize());
      TEST_EXIT_DBG(static_cast<int>( a.getSize()) >= a.getRowFESpace()->getAdmin()->getUsedSize())
336
	("a.size = %d too small: admin->sizeUsed = %d\n",
337
	a.getSize(), a.getRowFESpace()->getAdmin()->getUsedSize());
338
339
    }
    else  if (transpose == Transpose) {
340
      TEST_EXIT_DBG(static_cast<int>(x.getSize()) >= a.getRowFESpace()->getAdmin()->getUsedSize())
341
	("x.size = %d too small: admin->sizeUsed = %d\n",
342
343
	 x.getSize(), a.getRowFESpace()->getAdmin()->getUsedSize());
      TEST_EXIT_DBG(static_cast<int>(y.getSize()) >= a.getColFESpace()->getAdmin()->getUsedSize())
344
	("y.size = %d too small: admin->sizeUsed = %d\n",
345
346
	 y.getSize(), a.getColFESpace()->getAdmin()->getUsedSize());
      TEST_EXIT_DBG(static_cast<int>( a.getSize()) >= a.getRowFESpace()->getAdmin()->getUsedSize())
347
	("a.size = %d too small: admin->sizeUsed = %d\n",
348
	a.getSize(), a.getRowFESpace()->getAdmin()->getUsedSize());
349
350
351
352
353
    }

    ysize = y.getSize();

    typename DOFVector<T>::Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(&y), FREE_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
354
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator) { 
355
      *vecIterator = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
356
    }
357
358
359
360

    if (transpose == NoTranspose) {
      typename DOFVector<T>::Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(&y), USED_DOFS); 
      DOFMatrix::Iterator rowIterator(const_cast<DOFMatrix*>(&a), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
361
362
      for (vecIterator.reset(), rowIterator.reset(); !rowIterator.end();
	   ++rowIterator, ++vecIterator) { 
363
364
	sum = 0;
	row = &(a[rowIterator.getDOFIndex()]);
Thomas Witkowski's avatar
Thomas Witkowski committed
365
366
367
	for (std::vector<MatEntry>::iterator colIterator = rowIterator->begin();
	     colIterator != rowIterator->end();
	     colIterator++) {
368
369
370
371
372
373
374
375
376
377
	  jcol = colIterator->col;
	  if (jcol >= 0) { // entry used? 
	    sum += (static_cast<T>(colIterator->entry)) * x[jcol];
	  } else {
	    if (jcol == DOFMatrix::NO_MORE_ENTRIES)
	      break;
	  }
	}
	*vecIterator *= beta;
	*vecIterator += alpha * sum;
Thomas Witkowski's avatar
Thomas Witkowski committed
378
      }
379
380
    } else if (transpose == Transpose) {
      typename DOFVector<T>::Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(&y), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
381
      for (vecIterator.reset(); !vecIterator.end(); ++vecIterator) {
382
	*vecIterator  *= beta ; 
Thomas Witkowski's avatar
Thomas Witkowski committed
383
      }
384
385
386
    
      typename DOFVector<T>::Iterator xIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&x)), USED_DOFS);
      DOFMatrix::Iterator rowIterator(const_cast<DOFMatrix*>(&a), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
387
388
      for (xIterator.reset(), rowIterator.reset(); !rowIterator.end();
	   ++rowIterator, ++xIterator) { 
389
390
	ax = alpha * (*xIterator);
	row = &(a[rowIterator.getDOFIndex()]);
Thomas Witkowski's avatar
Thomas Witkowski committed
391
	for (std::vector<MatEntry>::iterator colIterator = rowIterator->begin();
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
	    colIterator != rowIterator->end();
	    colIterator++) {
	  jcol = colIterator->col;
	  if (jcol >= 0) // entry used?
	    y[jcol] += ax * (static_cast<T>(colIterator->entry));
	  else 
	    if (jcol == DOFMatrix::NO_MORE_ENTRIES) break;
	}
      }
    }
    else {
      ERROR_EXIT("transpose=%d\n", transpose);
    }
  }

  template<typename T>
  void DOFVector<T>::print() const
  {
410
    FUNCNAME("DOFVector<T>::print()");
411
    const DOFAdmin *admin = NULL;
412
    const char *format;
413

414
415
    if (feSpace) 
      admin = feSpace->getAdmin();
416
417

    MSG("Vec `%s':\n", name.c_str());
418
    int j = 0;
419
420
421
422
423
424
425
426
427
    if (admin) {
      if (admin->getUsedSize() > 100)
	format = "%s(%3d,%10.5le)";
      else if (admin->getUsedSize() > 10)
	format = "%s(%2d,%10.5le)";
      else
	format = "%s(%1d,%10.5le)";

      Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(this)), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
428
      for (vecIterator.reset(); !vecIterator.end(); ++vecIterator) {
429
	if ((j % 3) == 0) {
430
431
	  if (j) 
	    Msg::print("\n");
432
	  MSG(format, "", vecIterator.getDOFIndex(), *vecIterator);
Thomas Witkowski's avatar
Thomas Witkowski committed
433
	} else {
434
	  Msg::print(format, " ", vecIterator.getDOFIndex(), *vecIterator);
Thomas Witkowski's avatar
Thomas Witkowski committed
435
	}
436
	j++;
437
      }
438
      Msg::print("\n");
439
    } else {
440
441
	MSG("no DOFAdmin, print whole vector.\n");
    
442
443
444
445
446
447
	for (int i = 0; i < static_cast<int>( vec.size()); i++) {
	  if ((j % 3) == 0) {
	    if (j) 
	      Msg::print("\n");
	    MSG("(%d,%10.5e)",i,vec[i]);
	  } else {
448
	    Msg::print(" (%d,%10.5e)",i,vec[i]);
449
	  }
450
451
452
453
454
455
456
	  j++;
	}
	Msg::print("\n");
      }
    return;
  }

457
458
459
460
461
462
463
464
465
  template<typename T>
  int DOFVector<T>::calcMemoryUsage() const
  {
    int result = 0;
    result += sizeof(DOFVector<T>);
    result += vec.size() * sizeof(T);

    return result;
  }
466
467
468
469
470
471
472
473
474

  template<typename T>
  T DOFVectorBase<T>::evalUh(const DimVec<double>& lambda, 
			     DegreeOfFreedom* dof_indices)
  {
    BasisFunction* phi = const_cast<BasisFunction*>(this->getFESpace()->getBasisFcts());
    int numberOfBasFcts = phi->getNumber();
    T val = 0.0;

475
    for (int i = 0; i < numberOfBasFcts; i++)
476
477
478
479
480
481
482
483
      val += (*this)[dof_indices[i]]*(*phi->getPhi(i))(lambda);

    return val;
  }

  template<typename T>
  void DOFVector<T>::interpol(AbstractFunction<T, WorldVector<double> > *fct)
  {
484
    FUNCNAME("DOFVector<T>::interpol()");
485

486
487
488
    TEST_EXIT_DBG(fct)("no function to interpolate\n");

    interFct = fct;
489

490
491
492
493
494
    if (!this->getFESpace()) {
      MSG("no dof admin in vec %s, skipping interpolation\n",
	  this->getName().c_str());
      return;
    }
495

496
497
498
499
500
    if (!(this->getFESpace()->getAdmin())) {
      MSG("no dof admin in feSpace %s, skipping interpolation\n", 
	  this->getFESpace()->getName().c_str());
      return;
    }
501
  
502
503
504
505
506
    if (!(this->getFESpace()->getBasisFcts())) {
      MSG("no basis functions in admin of vec %s, skipping interpolation\n",
	  this->getName().c_str());
      return;
    }
507

508
509
510
511
512
    if (!(fct)) {
      MSG("function that should be interpolated only pointer to NULL, ");
      Msg::print("skipping interpolation\n");
      return;
    }
513
514
515

    traverseVector = this;

Thomas Witkowski's avatar
Thomas Witkowski committed
516
517
518
519
520
521
522
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(this->getFESpace()->getMesh(), -1,
					 Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
    while (elInfo) {
      interpolFct(elInfo);
      elInfo = stack.traverseNext(elInfo);
    }
523
524
525
526
527
528
529
  }

  template<typename T>
  int DOFVector<T>::interpolFct(ElInfo* elinfo)
  {
    const BasisFunction *basFct = traverseVector->getFESpace()->getBasisFcts();
    const DOFAdmin* admin = traverseVector->getFESpace()->getAdmin();
Thomas Witkowski's avatar
Thomas Witkowski committed
530
531
532
533
534
535
536
    const DegreeOfFreedom *dof =  basFct->getLocalIndices(const_cast<Element *>(elinfo->getElement()),
							  admin, NULL);
    const T *inter_val = const_cast<BasisFunction*>(basFct)->interpol(elinfo, 
								      0, 
								      NULL,
								      traverseVector->interFct, 
								      NULL);
537
538

    int number = basFct->getNumber();
Thomas Witkowski's avatar
Thomas Witkowski committed
539
    for (int i = 0; i < number; i++)
540
541
542
543
544
545
546
547
      (*traverseVector)[dof[i]] = inter_val[i];

    return 0;
  }

  template<typename T>
  double DOFVector<T>::Int(Quadrature* q) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
548
    FUNCNAME("DOFVector::Int()");
549
550
551
  
    Mesh* mesh = feSpace->getMesh();

Thomas Witkowski's avatar
Thomas Witkowski committed
552
553
554
555
    if (!q) {
      int deg = 2 * feSpace->getBasisFcts()->getDegree();
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577

    quad_fast = FastQuadrature::provideFastQuadrature(feSpace->getBasisFcts(), *q, INIT_PHI);
    norm = 0.0;
    traverseVector = const_cast<DOFVector<T>*>(this);

    mesh->traverse(-1, 
		   Mesh::CALL_LEAF_EL|
		   Mesh::FILL_COORDS |
		   Mesh::FILL_DET, 
		   Int_fct);

    return norm;  
  }

  template<typename T>
  int DOFVector<T>::Int_fct(ElInfo *elinfo)
  {
    double det, normT;
    const T *uh_vec;
    int iq;

    det = elinfo->getDet();
578
    uh_vec = traverseVector->getVecAtQPs(elinfo, NULL, quad_fast, NULL);
579
580

    int numPoints = quad_fast->getNumPoints();
581
582
583
    for (normT = iq = 0; iq < numPoints; iq++) {
      normT += quad_fast->getWeight(iq)*(uh_vec[iq]);
    }
584
585
586
587
588
589
590
591
    norm += det*normT;

    return 0;
  }

  template<typename T>
  double DOFVector<T>::L1Norm(Quadrature* q) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
592
    FUNCNAME("DOFVector::L1Norm()");
593
  
Thomas Witkowski's avatar
Thomas Witkowski committed
594
595
596
597
    if (!q) {
      int deg = 2 * feSpace->getBasisFcts()->getDegree();
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
598
599
600
601

    quad_fast = FastQuadrature::provideFastQuadrature(feSpace->getBasisFcts(),*q,INIT_PHI);
    norm = 0.0;
    traverseVector = const_cast<DOFVector<T>*>(this);
Thomas Witkowski's avatar
Thomas Witkowski committed
602
    Mesh* mesh = feSpace->getMesh();
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620

    mesh->traverse(-1, 
		   Mesh::CALL_LEAF_EL|
		   Mesh::FILL_COORDS |
		   Mesh::FILL_DET, 
		   L1Norm_fct);

    return norm;  
  }

  template<typename T>
  int DOFVector<T>::L1Norm_fct(ElInfo *elinfo)
  {
    double det, normT;
    const T *uh_loc, *uh_vec;
    int iq;

    det = elinfo->getDet();
621
    uh_vec = traverseVector->getVecAtQPs(elinfo, NULL, quad_fast, NULL);
622
623

    int numPoints = quad_fast->getNumPoints();
624
625
626
    for (normT = iq = 0; iq < numPoints; iq++) {
      normT += quad_fast->getWeight(iq)*abs(uh_vec[iq]);
    }
627
628
629
630
631
632
633
634
635
    norm += det*normT;

    return 0;
  }


  template<typename T>
  double DOFVector<T>::L2NormSquare(Quadrature* q) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
636
    FUNCNAME("DOFVector::L2NormSquare()");
637

Thomas Witkowski's avatar
Thomas Witkowski committed
638
639
640
641
    if (!q) {
      int deg = 2 * feSpace->getBasisFcts()->getDegree();
      q = Quadrature::provideQuadrature(this->dim, deg);
    }
642
643
644
645
646
647
648

    quad_fast = FastQuadrature::provideFastQuadrature(feSpace->getBasisFcts(), 
						      *q, 
						      INIT_PHI);

    norm = 0.0;
    traverseVector = const_cast<DOFVector<T>*>(this);
Thomas Witkowski's avatar
Thomas Witkowski committed
649
650
    Mesh* mesh = feSpace->getMesh();
    mesh->traverse(-1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | Mesh::FILL_DET, L2Norm_fct);
651
652
653
654
655
656
657
658
659
660
661
662

    return norm;  
  }

  template<typename T>
  int DOFVector<T>::L2Norm_fct(ElInfo *elinfo)
  {
    double det, normT;
    const T *uh_vec;
    int iq;

    det = elinfo->getDet();
663
    uh_vec = traverseVector->getVecAtQPs(elinfo, NULL, quad_fast, NULL);
664
665

    int numPoints = quad_fast->getNumPoints();
666
667
668
    for (normT = iq = 0; iq < numPoints; iq++) {
      normT += quad_fast->getWeight(iq)*sqr(uh_vec[iq]);
    }
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
    norm += det*normT;

    return 0;
  }

  template<typename T>
  int DOFVector<T>::H1Norm_fct(ElInfo *elinfo)
  {
    double norm2, normT;
    const WorldVector<T> *grduh_vec;
    int iq, j;

    double det = elinfo->getDet();

    grduh_vec = traverseVector->getGrdAtQPs(elinfo, NULL, quad_fast, NULL);

    int dimOfWorld = Global::getGeo(WORLD);

    int numPoints = quad_fast->getNumPoints();
    for (normT = iq = 0; iq < numPoints; iq++)
      {
	for (norm2 = j = 0; j < dimOfWorld; j++)
	  norm2 += sqr(grduh_vec[iq][j]);

	normT += quad_fast->getWeight(iq)*norm2;
      }
    norm += det*normT;

    return 0;
  }


  template<typename T>
  double DOFVector<T>::H1NormSquare(Quadrature *q) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
704
    FUNCNAME("DOFVector::H1NormSquare()");
705

Thomas Witkowski's avatar
Thomas Witkowski committed
706
707
708
709
710
711
712
    if (!q) {
      int deg = 2 * feSpace->getBasisFcts()->getDegree() - 2;
      q = Quadrature::provideQuadrature(this->dim, deg);
    }

    quad_fast = FastQuadrature::provideFastQuadrature(feSpace->getBasisFcts(), 
						      *q, INIT_GRD_PHI);
713
714

    norm = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
715
716
    traverseVector = const_cast<DOFVector<T>*>(this); 
    Mesh *mesh = feSpace->getMesh();
717
718
719
720
721
722
723
724
725
726

    mesh->traverse(-1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | 
		   Mesh::FILL_DET | Mesh::FILL_GRD_LAMBDA,
		   H1Norm_fct);

    return norm;
  }

  template<typename T>
  void DOFVector<T>::compressDOFIndexed(int first, int last, 
727
					std::vector<DegreeOfFreedom> &newDOF)
728
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
729
730
    for (int i = first; i <= last; i++)
      if (newDOF[i] >= 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
731
	vec[newDOF[i]] = vec[i];
732
733
734
735
736
737
  }

  template<typename T>
  Flag DOFVectorBase<T>::getAssembleFlag()
  {
    Flag fillFlag(0);
Thomas Witkowski's avatar
Thomas Witkowski committed
738
739
740
    
    for (std::vector<Operator*>::iterator op = this->operators.begin(); 
	 op != this->operators.end(); ++op) {
741
742
      fillFlag |= (*op)->getFillFlag();
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
743

744
745
746
    return fillFlag;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
747
748
749
750
751
752
753
754
755
756
  template<typename T>
  void DOFVectorBase<T>::finishAssembling()
  {
    // call the operatos cleanup procedures
    for (std::vector<Operator*>::iterator it = this->operators.begin();
	 it != this->operators.end(); ++it) {     
      (*it)->finishAssembling();
    }
  }

757
  template<typename T>
758
759
  DOFVector<T>& DOFVector<T>::operator=(const DOFVector<T>& rhs )
  {
760
    feSpace = rhs.feSpace;
761
    DOFVectorBase<T>::feSpace = rhs.feSpace;
762
    this->nBasFcts = rhs.nBasFcts;
763
    vec = rhs.vec;
764
    this->elementVector = new ElementVector(this->nBasFcts);
765
766
767
768
    interFct = rhs.interFct;
    coarsenOperation = rhs.coarsenOperation;
    this->operators = rhs.operators;
    this->operatorFactor = rhs.operatorFactor;
769
    if (rhs.boundaryManager) {
770
771
      if (this->boundaryManager) 
	delete this->boundaryManager; 
772
      this->boundaryManager = new BoundaryManager(*rhs.boundaryManager);
773

774
      //    boundaryManager->setDOFVector(this);
775
    } else {
776
      this->boundaryManager = NULL;
777
778
    }

779
780
781
782
783
784
    return *this;
  }

  template<typename T>
  const DOFVector<T>& operator*=(DOFVector<T>& x, T scal)
  {
785
    FUNCNAME("DOFVector<T>::operator*=(DOFVector<T>& x, T scal)");
786

787
788
    TEST_EXIT_DBG(x.getFESpace() && x.getFESpace()->getAdmin())
      ("pointer is NULL: %8X, %8X\n", x.getFESpace(), x.getFESpace()->getAdmin());
789
790

    typename DOFVector<T>::Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(&x), USED_DOFS);
791
    for (vecIterator.reset(); !vecIterator.end(); ++vecIterator) {
792
      (*vecIterator) *= scal; 
Thomas Witkowski's avatar
Thomas Witkowski committed
793
    }
794

795
796
797
798
799
800
801
    return x;
  }


  template<typename T>
  const DOFVector<T>& operator+=(DOFVector<T>& x, const DOFVector<T>& y)
  {
802
803
804
    FUNCNAME("DOFVector<T>::operator+=(DOFVector<T>& x, const DOFVector<T>& y)");
    
    TEST_EXIT_DBG(x.getFESpace() && y.getFESpace())
805
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());
806
807
    TEST_EXIT_DBG(x.getFESpace()->getAdmin() &&
	      (x.getFESpace()->getAdmin() == y.getFESpace()->getAdmin()))
808
809
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
810
811
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");
    
812
813
    typename DOFVector<T>::Iterator xIterator(dynamic_cast<DOFIndexed<T>*>(&x), USED_DOFS);
    typename DOFVector<T>::Iterator yIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&y)), USED_DOFS);
814
815
816
817
818
819
    for (xIterator.reset(), yIterator.reset();
	 !xIterator.end();
	 ++xIterator, ++yIterator) {     
      *xIterator += *yIterator; 
    }

820
821
822
823
824
825
    return x;
  }

  template<typename T>
  const DOFVector<T>& operator-=(DOFVector<T>& x, const DOFVector<T>& y)
  {
826
    FUNCNAME("DOFVector<T>::operator-=(DOFVector<T>& x, const DOFVector<T>& y)");
827

828
    TEST_EXIT_DBG(x.getFESpace() && y.getFESpace())
829
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());
830
831
    TEST_EXIT_DBG(x.getFESpace()->getAdmin() &&
	      (x.getFESpace()->getAdmin() == y.getFESpace()->getAdmin()))
832
833
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
834
835
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");
    
836
837
    typename DOFVector<T>::Iterator xIterator(dynamic_cast<DOFIndexed<T>*>(&x), USED_DOFS);
    typename DOFVector<T>::Iterator yIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&y)), USED_DOFS);
838
839
840
841
842
    for (xIterator.reset(), yIterator.reset();
	 !xIterator.end();
	 ++xIterator, ++yIterator) {
      *xIterator -= *yIterator; 
    }
843
844
845
846
847
848
    return x;
  }

  template<typename T>
  const DOFVector<T>& operator*=(DOFVector<T>& x, const DOFVector<T>& y)
  {
849
850
851
    FUNCNAME("DOFVector<T>::operator*=(DOFVector<T>& x, const DOFVector<T>& y)");
    
    TEST_EXIT_DBG(x.getFESpace() && y.getFESpace())
852
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());
853
854
    TEST_EXIT_DBG(x.getFESpace()->getAdmin() &&
	      (x.getFESpace()->getAdmin() == y.getFESpace()->getAdmin()))
855
856
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
857
858
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");
    
859
860
    typename DOFVector<T>::Iterator xIterator(dynamic_cast<DOFIndexed<T>*>(&x), USED_DOFS);
    typename DOFVector<T>::Iterator yIterator(dynamic_cast<DOFIndexed<T>*>(const_cast<DOFVector<T>*>(&y)), USED_DOFS);
861
862
863
864
865
    for (xIterator.reset(), yIterator.reset();
	 !xIterator.end();
	 ++xIterator, ++yIterator) {
      *xIterator *= *yIterator; 
    }
866
867
868
869
870
871
872
873
874

    return x;
  }

  template<typename T>
  T operator*(DOFVector<T>& x, DOFVector<T>& y)
  {
    FUNCNAME("DOFVector<T>::operator*(DOFVector<T>& x, DOFVector<T>& y)");

875
    TEST_EXIT_DBG(x.getFESpace() && y.getFESpace())
876
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());
877
878
    TEST_EXIT_DBG(x.getFESpace()->getAdmin() && 
	      (x.getFESpace()->getAdmin() == y.getFESpace()->getAdmin()))
879
880
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
881
    TEST_EXIT_DBG(x.getSize() == y.getSize())("different sizes\n");
882

883
    T dot = 0;    
884
885
    typename DOFVector<T>::Iterator xIterator(dynamic_cast<DOFIndexed<T>*>(&x), USED_DOFS);
    typename DOFVector<T>::Iterator yIterator(dynamic_cast<DOFIndexed<T>*>(&y), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
886
    for (xIterator.reset(), yIterator.reset(); !xIterator.end(); 
887
888
	 ++xIterator, ++yIterator) {
      dot += (*xIterator) * (*yIterator);
Thomas Witkowski's avatar
Thomas Witkowski committed
889
    }
890
891
892
893
894
895
896
897
898
899
900
901
902

    return(dot);
  }

  template<typename T>
  void mv(MatrixTranspose transpose,
	  const DOFMatrix &a, 
	  const DOFVector<T>&x,
	  DOFVector<T> &result,
	  bool add)
  {
    FUNCNAME("DOFVector<T>::mv()");

903
    TEST_EXIT_DBG(a.getRowFESpace() && a.getColFESpace() && x.getFESpace() && result.getFESpace())
904
905
      ("getFESpace() is NULL: %8X, %8X, %8X, %8X\n", 
       a.getRowFESpace(), a.getColFESpace(), x.getFESpace(), result.getFESpace());
906
907
908
909
910
911
912
913
    
    TEST_EXIT_DBG((a.getRowFESpace()->getAdmin() && a.getColFESpace()->getAdmin()) && 
	      (((transpose == NoTranspose) && 
		(a.getColFESpace()->getAdmin() == x.getFESpace()->getAdmin()) && 
		(a.getRowFESpace()->getAdmin() == result.getFESpace()->getAdmin())) ||
	       ((transpose == Transpose) && 
		(a.getRowFESpace()->getAdmin() == x.getFESpace()->getAdmin()) && 
		(a.getColFESpace()->getAdmin() == result.getFESpace()->getAdmin()))))
914
915
916
      ("no admin or different admins: %8X, %8X, %8X, %8X\n",
       a.getRowFESpace()->getAdmin(), a.getColFESpace()->getAdmin(), 
       x.getFESpace()->getAdmin(), result.getFESpace()->getAdmin());
917
    
918
    if (transpose == NoTranspose) {
919
      TEST_EXIT_DBG(static_cast<int>(x.getSize()) >= a.getColFESpace()->getAdmin()->getUsedSize())
920
	("x.size = %d too small: admin->sizeUsed = %d\n",
921
922
	 x.getSize(), a.getColFESpace()->getAdmin()->getUsedSize());
      TEST_EXIT_DBG(static_cast<int>( result.getSize()) >= a.getRowFESpace()->getAdmin()->getUsedSize())
923
	("size = %d too small: admin->sizeUsed = %d\n",
924
925
	 result.getSize(), a.getRowFESpace()->getAdmin()->getUsedSize());
      TEST_EXIT_DBG(static_cast<int>( a.getSize()) >= a.getRowFESpace()->getAdmin()->getUsedSize())
926
	("a.size = %d too small: admin->sizeUsed = %d\n",
927
928
	 a.getSize(), a.getRowFESpace()->getAdmin()->getUsedSize());
      
929

Thomas Witkowski's avatar
Thomas Witkowski committed
930
931
      typename DOFVector<T>::Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(&result), USED_DOFS); 
      DOFMatrix::Iterator rowIterator(const_cast<DOFMatrix*>(&a), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
932
      for (vecIterator.reset(), rowIterator.reset(); !rowIterator.end();
Thomas Witkowski's avatar
Thomas Witkowski committed
933
934
935
936
937
      	  ++rowIterator, ++vecIterator) { 
      	
	double sum = 0;
      	if (!add) 
	  *vecIterator = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
938
939
940
      	for (std::vector<MatEntry>::iterator colIterator = rowIterator->begin();
	     colIterator != rowIterator->end();
	     colIterator++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
941
942
943
944
945
946
947
948
949
950
      	  
	  int jcol = colIterator->col;
      	  if (jcol >= 0) { // entry used? 
      	    sum += (static_cast<double>(colIterator->entry)) * x[jcol];
      	  } else {
      	    if (jcol == DOFMatrix::NO_MORE_ENTRIES)
      	      break;
      	  }
      	}
      	*vecIterator += sum;
951
952
953
      }

    } else if (transpose == Transpose) {
954
      TEST_EXIT_DBG(static_cast<int>(x.getSize()) >= a.getRowFESpace()->getAdmin()->getUsedSize())
955
	("x.size = %d too small: admin->sizeUsed = %d\n",
956
957
	 x.getSize(), a.getRowFESpace()->getAdmin()->getUsedSize());
      TEST_EXIT_DBG(static_cast<int>( result.getSize()) >= a.getColFESpace()->getAdmin()->getUsedSize())
958
	("size = %d too small: admin->sizeUsed = %d\n",
959
960
	 result.getSize(), a.getColFESpace()->getAdmin()->getUsedSize());
      TEST_EXIT_DBG(static_cast<int>( a.getSize()) >= a.getRowFESpace()->getAdmin()->getUsedSize())
961
	("a.size = %d too small: admin->sizeUsed = %d\n",
962
963
	 a.getSize(), a.getRowFESpace()->getAdmin()->getUsedSize());
      
Thomas Witkowski's avatar
Thomas Witkowski committed
964
      if (!add) {
965
	typename DOFVector<T>::Iterator vecIterator(dynamic_cast<DOFIndexed<T>*>(&result), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
966
	for (vecIterator.reset(); !vecIterator.end(); ++vecIterator) {
967
	  *vecIterator = 0; 
Thomas Witkowski's avatar
Thomas Witkowski committed
968
	}
969
970
971
972
      }

      typename DOFVector<T>::Iterator xIterator(const_cast<DOFIndexed<T>*>(dynamic_cast<const DOFIndexed<T>*>(&x)), USED_DOFS);
      DOFMatrix::Iterator rowIterator(const_cast<DOFMatrix*>(&a), USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
973
974
      for (xIterator.reset(), rowIterator.reset(); !rowIterator.end();
	   ++rowIterator, ++xIterator) { 
975
	T ax = (*xIterator);
976
	for (std::vector<MatEntry>::iterator colIterator = rowIterator->begin();
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
	    colIterator != rowIterator->end();
	    colIterator++) {
	  int jcol = colIterator->col;
	  if (jcol >= 0) // entry used?
	    result[jcol] += (static_cast<double>(colIterator->entry)) * ax;
	  else 
	    if (jcol == DOFMatrix::NO_MORE_ENTRIES) break;
	}
      }
    } else {
      ERROR_EXIT("transpose = %d\n", transpose);
    }
  }

  template<typename T>
  void axpy(double alpha, 
	    const DOFVector<T>& x, 
	    DOFVector<T>& y)
  {
    FUNCNAME("DOFVector<T>::axpy()");

998
    TEST_EXIT_DBG(x.getFESpace() && y.getFESpace())
999
1000
      ("feSpace is NULL: %8X, %8X\n", x.getFESpace(), y.getFESpace());

1001
1002
    TEST_EXIT_DBG(x.getFESpace()->getAdmin() && 
	      (x.getFESpace()->getAdmin() == y.getFESpace()->getAdmin()))
1003
1004
      ("no admin or different admins: %8X, %8X\n",
       x.getFESpace()->getAdmin(), y.getFESpace()->getAdmin());
1005
    TEST_EXIT_DBG(static_cast<int>(x.getSize()) >= x.getFESpace()->getAdmin()->getUsedSize())
1006
      ("size = %d too small: admin->size = %d\n", x.getSize(), 
1007
1008
       x.getFESpace()->getAdmin()->getUsedSize());
    TEST_EXIT_DBG(static_cast<int>(y.getSize()) >= x.getFESpace()->getAdmin()->