DOFMatrix.cc 24.3 KB
Newer Older
1 2
#include <algorithm>
#include <png.h>
3 4

#include "DOFMatrix.h"
5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30
#include "QPsiPhi.h"
#include "BasisFunction.h"
#include "Boundary.h"
#include "DOFAdmin.h"
#include "ElInfo.h"
#include "FiniteElemSpace.h"
#include "Mesh.h"
#include "DOFVector.h"
#include "Operator.h"
#include "BoundaryCondition.h"
#include "BoundaryManager.h"
#include "ElementMatrix.h"
#include "Assembler.h"

namespace AMDiS {

  DOFMatrix *DOFMatrix::traversePtr = NULL;

  DOFMatrix::DOFMatrix()
  {
    rowFESpace = NULL;
    colFESpace = NULL;
  }

  DOFMatrix::DOFMatrix(const FiniteElemSpace* rowFESpace_,
		       const FiniteElemSpace* colFESpace_,
31
		       std::string name_)
32 33 34 35 36 37 38
    : rowFESpace(rowFESpace_),
      colFESpace(colFESpace_),
      name(name_), 
      coupleMatrix(false)
  {
    TEST_EXIT(rowFESpace)("no rowFESpace\n");
  
39
    if (!colFESpace) {
40 41 42 43
      colFESpace = rowFESpace;
    }

    if (rowFESpace && rowFESpace->getAdmin())
Thomas Witkowski's avatar
Thomas Witkowski committed
44
      (const_cast<DOFAdmin*>(rowFESpace->getAdmin()))->addDOFIndexed(this);
45

Thomas Witkowski's avatar
Thomas Witkowski committed
46
    boundaryManager = NEW BoundaryManager(rowFESpace_);
Thomas Witkowski's avatar
Thomas Witkowski committed
47 48
    elementMatrix = NEW ElementMatrix(rowFESpace->getBasisFcts()->getNumber(),
				      colFESpace->getBasisFcts()->getNumber());
Thomas Witkowski's avatar
Thomas Witkowski committed
49 50

    applyDBCs.clear();
51 52 53 54

#ifdef HAVE_PARALLEL_AMDIS
    applicationOrdering = NULL;
#endif
55 56 57
  }

  DOFMatrix::DOFMatrix(const DOFMatrix& rhs)
Thomas Witkowski's avatar
Thomas Witkowski committed
58
    : name(rhs.name + "copy")
59
  {
60
    *this = rhs;
61 62 63 64 65 66
    if (rowFESpace && rowFESpace->getAdmin())
      (const_cast<DOFAdmin*>( rowFESpace->getAdmin()))->addDOFIndexed(this);
  }

  DOFMatrix::~DOFMatrix()
  {
67 68
    FUNCNAME("DOFMatrix::~DOFMatrix()");

69 70
    matrix.clear();

71
    if (rowFESpace && rowFESpace->getAdmin()) {
72 73
      (const_cast<DOFAdmin*>(rowFESpace->getAdmin()))->removeDOFIndexed(this);
    }  
Thomas Witkowski's avatar
Thomas Witkowski committed
74

75 76 77
    if (boundaryManager)
      DELETE boundaryManager;

Thomas Witkowski's avatar
Thomas Witkowski committed
78 79
    if (elementMatrix)
      DELETE elementMatrix;
80 81 82 83 84 85 86 87 88 89 90 91 92
  }

  void DOFMatrix::print() const
  {
    FUNCNAME("DOFMatrix::print()");

    int sizeUsed = rowFESpace->getAdmin()->getUsedSize();

    if (static_cast<int>(matrix.size()) < sizeUsed) {
      WARNING("DOFMatrix not yet initialized\n");
      return;
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
93 94
    for (int i = 0; i < sizeUsed; i++) { 
      DOFMatrix::MatrixRow row = matrix[i];
95 96
      MSG("row %3d:",i);
      int rowSize = static_cast<int>( row.size());
Thomas Witkowski's avatar
Thomas Witkowski committed
97 98 99
      for (int j = 0; j < rowSize; j++) {
	int jcol = row[j].col;
	if (entryUsed(i, j)) {
100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135
	  Msg::print(" (%3d,%20.17lf)", jcol, row[j].entry);
	}
      }
      Msg::print("\n");
    }
  }

  void DOFMatrix::printRow(int i) const
  {
    FUNCNAME("DOFMatrix::printRow()");

    int sizeUsed = rowFESpace->getAdmin()->getUsedSize();

    if (static_cast<int>(matrix.size()) < sizeUsed) {
      WARNING("DOFMatrix not yet initialized\n");
      return;
    }

    DOFMatrix::MatrixRow row = matrix[i];
    MSG("row %3d:",i);
    int rowSize = static_cast<int>( row.size());
    for (int j = 0; j < rowSize; j++) {
      int jcol = row[j].col;
      if (entryUsed(i,j)) {
	Msg::print(" (%3d,%20.17lf)", jcol, row[j].entry);
      }
    }
    Msg::print("\n");   
  }

  /****************************************************************************/
  /*  clear: remove all entries from dof_matrix                               */
  /****************************************************************************/

  void DOFMatrix::clear()
  {
136 137
    int mSize = static_cast<int>(matrix.size());
    for (int i = 0; i < mSize; i++) {
138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229
      matrix[i].resize(0);
    }
  }

  bool DOFMatrix::symmetric()
  {
    FUNCNAME("DOFMatrix::symmetric()");

    DegreeOfFreedom row, col;
    double entry, tol = 1e-5;

    DOFMatrix::Iterator matrixRow(this, USED_DOFS);

    for (matrixRow.reset(); !matrixRow.end(); ++matrixRow) {
      row = matrixRow.getDOFIndex();
      int rowSize = matrixRow->size();
      for (int i = 0; i < rowSize; i++) {
	col = (*matrixRow)[i].col;
	entry = (*matrixRow)[i].entry;
	if (abs(entry - logAcc(col, row)) > tol) {
	  MSG("matrix[%d][%d] = %e, matrix[%d][%d] = %e\n",
	      row, col, entry, col, row, logAcc(col, row));

	  return false;
	}
      }
    }

    return true;
  }

  void DOFMatrix::test()
  {
    FUNCNAME("DOFMatrix::test()");

    int non_symmetric = 0, found = 0;

    /* test symmetry */
    for (int i = 0; i < static_cast<int>(matrix.size()); i++) { 
      double sum = 0.0;
      DOFMatrix::MatrixRow *row = &matrix[i];
      for (int j = 0; j < static_cast<int>(row->size()); j++) {
	int jcol = (*row)[j].col;
	if (entryUsed(i, j)) {
	  found = 0;
	  if ((*row)[j].entry != (*row)[j].entry) {
	    MSG("mat[%d,%d]=%10.5e ???\n", i, jcol, (*row)[j].entry);
	    WAIT;
	  }
	  DOFMatrix::MatrixRow *row2 = &matrix[jcol];
	  for (int k = 0; k < static_cast<int>(row->size()); k++) {
	    int kcol = (*row)[k].col;
	    if (entryUsed(jcol, k)) {
	      if (kcol == i) {
		found = 1;
		if (abs((*row2)[k].entry - (*row)[j].entry) > 1.E-5) {
		  non_symmetric = 1;
		  MSG("mat[%d,%d]=%10.5e != mat[%d,%d]=%10.5e\n",
		      i, jcol, (*row)[j].entry, jcol, i, (*row2)[k].entry);
		}
		row2 = NULL;
		break;
	      }
	    }
	  }
	  if (!found) {
	    non_symmetric = 1;
	    MSG("mat[%d,%d] not found\n", jcol, i);
	  }
	}
      }
      if (abs(sum) > 1.E-5) {
	MSG("Zeilensumme[%d] = %10.5e\n", i, sum);
      }
    }

    if (non_symmetric) {
      MSG("matrix `%s' not symmetric.\n", name.data());
    } else {
      MSG("matrix `%s' is symmetric.\n", name.data());
    }
  }


  DOFMatrix& DOFMatrix::operator=(const DOFMatrix& rhs)
  {
    rowFESpace = rhs.rowFESpace;
    colFESpace = rhs.colFESpace;
    operators = rhs.operators;
    operatorFactor = rhs.operatorFactor;
    matrix = rhs.matrix;
    if (rhs.boundaryManager) {
Thomas Witkowski's avatar
Thomas Witkowski committed
230
      boundaryManager = NEW BoundaryManager(*rhs.boundaryManager);
231
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
232 233 234 235 236 237 238
      boundaryManager = NULL;
    }
    if (rhs.elementMatrix) {
      elementMatrix = NEW ElementMatrix(rowFESpace->getBasisFcts()->getNumber(),
					colFESpace->getBasisFcts()->getNumber());
    } else {
      elementMatrix = NULL;
239 240 241 242 243 244 245 246 247 248 249 250 251
    }

    return *this;
  }


  void DOFMatrix::addElementMatrix(double sign, 
				   const ElementMatrix &elMat, 
				   const BoundaryType *bound,
				   bool add)
  {
    FUNCNAME("DOFMatrix::addElementMatrix");

Thomas Witkowski's avatar
Thomas Witkowski committed
252
    DegreeOfFreedom row;
253

Thomas Witkowski's avatar
Thomas Witkowski committed
254 255
    int nRow = elMat.rowIndices.getSize();
    int nCol = elMat.colIndices.getSize();
256

Thomas Witkowski's avatar
Thomas Witkowski committed
257
    for (int i = 0; i < nRow; i++)  {   // for all rows of element matrix
258 259 260 261 262
      row = elMat.rowIndices[i];
      BoundaryCondition *condition = 
	bound ? boundaryManager->getBoundaryCondition(bound[i]) : NULL;

      if (condition && condition->isDirichlet()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
263
	applyDBCs.insert(static_cast<int>(row));
Thomas Witkowski's avatar
Thomas Witkowski committed
264
      } else {     
Thomas Witkowski's avatar
Thomas Witkowski committed
265 266 267
	for (int j = 0; j < nCol; j++) {  // for all columns
	  addSparseDOFEntry(sign, row, elMat.colIndices[j], 
			    elMat[i][j], add); 
268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289
	}
      }
    }
  }

  double DOFMatrix::logAcc(DegreeOfFreedom a,DegreeOfFreedom b) const
  {
    int j;

    for (j = 0; j < static_cast<int>(matrix[a].size()); j++) 
      if (b == matrix[a][j].col) 
	break;

    return (j == static_cast<int>(matrix[a].size())) ? 0.0 : matrix[a][j].entry;
  }

  void DOFMatrix::changeColOfEntry(DegreeOfFreedom a,
				   DegreeOfFreedom b,
				   DegreeOfFreedom c)
  {
    int j;

Thomas Witkowski's avatar
Thomas Witkowski committed
290 291 292 293 294 295
    for (j = 0; j<static_cast<int>(matrix[a].size());j++) 
      if (b == matrix[a][j].col) 
	break;

    if (j != static_cast<int>(matrix[a].size())) 
      matrix[a][j].col = c;
296 297
  }

298 299 300
  void DOFMatrix::addSparseDOFEntry(double sign, int irow, 
				    int jcol, double entry,
				    bool add)
301
  {
302
    FUNCNAME("DOFMatrix::addSparseDOFEntry()");
303 304 305

    MatrixRow *row = &(matrix[irow]);

306
    if (add && !entry) 
307
      return;
308

309
    int freeCol = -1;
Thomas Witkowski's avatar
Thomas Witkowski committed
310
    int rowSize = static_cast<int>(row->size());
311

312 313 314
    TEST_EXIT_DBG(jcol >= 0 && 
		  jcol < colFESpace->getAdmin()->getUsedSize())
      ("Column index %d out of range 0-%d\n", jcol, colFESpace->getAdmin()->getUsedSize() - 1);
315 316

    // first entry is diagonal entry
317
    if ((rowSize == 0) && (colFESpace == rowFESpace)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
318 319 320 321
      MatEntry newEntry = {irow, 0.0};
      row->push_back(newEntry);
      rowSize = 1;
    }
322

323
    int i;
324
    // search jcol
Thomas Witkowski's avatar
Thomas Witkowski committed
325
    for (i = 0; i < rowSize; i++) {
326
      // jcol found ?
Thomas Witkowski's avatar
Thomas Witkowski committed
327
      if ((*row)[i].col == jcol) {
328 329 330
	break;
      }
      // remember free entry
Thomas Witkowski's avatar
Thomas Witkowski committed
331
      if ((*row)[i].col == UNUSED_ENTRY) {
332 333 334
	freeCol = i;
      }
      // no more entries
Thomas Witkowski's avatar
Thomas Witkowski committed
335
      if ((*row)[i].col == NO_MORE_ENTRIES) {
336
	freeCol = i;
Thomas Witkowski's avatar
Thomas Witkowski committed
337
	if (rowSize > i+1) {
Thomas Witkowski's avatar
Thomas Witkowski committed
338
	  (*row)[i + 1].entry = NO_MORE_ENTRIES;
339 340 341 342 343 344
	}
	break;
      }
    }

    // jcol found?
Thomas Witkowski's avatar
Thomas Witkowski committed
345
    if (i < rowSize) {
Thomas Witkowski's avatar
Thomas Witkowski committed
346
      if (!add) {
Thomas Witkowski's avatar
Thomas Witkowski committed
347
	(*row)[i].entry = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
348 349
      }

350 351
      (*row)[i].entry += sign * entry;
    } else {
352
      if (freeCol == -1) {
353 354 355 356
	MatEntry newEntry = {jcol, sign * entry};
	row->push_back(newEntry);
      } else {
	(*row)[freeCol].col = jcol;
Thomas Witkowski's avatar
Thomas Witkowski committed
357 358
	if (!add) 
	  (*row)[freeCol].entry = 0.0;
359 360
	(*row)[freeCol].entry += sign * entry;
      }
361
    } 
362 363 364 365 366 367 368 369 370 371 372 373 374 375 376
  }

  void DOFMatrix::addMatEntry(int row, MatEntry entry)
  {
    matrix[row].push_back(entry);
  }

  void DOFMatrix::addMatEntry(int row,  DegreeOfFreedom col, double value)
  {
    MatEntry entry;
    entry.col = col;
    entry.entry = value;
    matrix[row].push_back(entry);
  }

377
  void DOFMatrix::addRow(std::vector<MatEntry> row)
378 379 380 381 382
  {
    matrix.push_back(row);
  }

  void DOFMatrix::compressDOFIndexed(int first, int last, 
383
				     std::vector<DegreeOfFreedom> &newDOF)
384
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
385 386 387
    for (int i = first; i <= last; i++) {
      if (newDOF[i] >= 0) {
	matrix[newDOF[i]].swap(matrix[i]);
388 389 390 391
	matrix[i].resize(0);
      }
    }
    int usedSize = rowFESpace->getAdmin()->getUsedSize();
Thomas Witkowski's avatar
Thomas Witkowski committed
392 393
    for (int i = 0; i < usedSize; i++) {
      std::vector<MatEntry> *row = reinterpret_cast<std::vector<MatEntry>*>(&(matrix[i])); 
394
      int rowSize = static_cast<int>(row->size());
Thomas Witkowski's avatar
Thomas Witkowski committed
395 396 397
      for (int j = 0; j < rowSize; j++) {
	int col = (*row)[j].col;
	if (entryUsed(i, j)) 
Thomas Witkowski's avatar
Thomas Witkowski committed
398
	  (*row)[j].col = newDOF[col];
399 400 401 402 403 404
      }    
    }
  }

  void DOFMatrix::freeDOFContent(int index)
  {
405
    int col = 0;
406 407

    if (0 < matrix[index].size()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
408 409
      int nRows = matrix.size();

410 411
      // for all columns in this row
      int size = static_cast<int>(matrix[index].size());
412
      for (int i = 0; i < size; i++) {
413
	// if entry is used
414
	if (entryUsed(index, i)) {
415 416
	  // get column of this entry
	  col = matrix[index][i].col;
Thomas Witkowski's avatar
Thomas Witkowski committed
417
	  if ((col != index) && (col < nRows)) {  // remove symmetric entry if exists
418
	    int colsize = static_cast<int>(matrix[col].size());
419 420
	    for (int j = 0; j < colsize; j++) {
	      int col2 = matrix[col][j].col;
421 422
	      if (col2 == index) {
		matrix[col][j].col = DOFMatrix::UNUSED_ENTRY;
423
	      } else if (col2 == DOFMatrix::NO_MORE_ENTRIES) {
424 425 426 427
		break;
	      }
	    }
	  }
428
	} else if (col == DOFMatrix::NO_MORE_ENTRIES) {
429 430 431 432 433 434 435 436
	  break;
	}
      }
      matrix[index].resize(0);	
    }
  }


437
  void DOFMatrix::assemble(double factor, ElInfo *elInfo, const BoundaryType *bound)
438
  {
439
    FUNCNAME("DOFMatrix::assemble()");
440

441 442
    if (operators.size() == 0)
	return;
443

Thomas Witkowski's avatar
Thomas Witkowski committed
444 445
    operators[0]->getAssembler(omp_get_thread_num())->
      initElementMatrix(elementMatrix, elInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
446

447 448 449
    std::vector<Operator*>::iterator it = operators.begin();
    std::vector<double*>::iterator factorIt = operatorFactor.begin();
    for (; it != operators.end(); ++it, ++factorIt) {
450
      if ((*it)->getNeedDualTraverse() == false) {
Thomas Witkowski's avatar
Thomas Witkowski committed
451
	  (*it)->getElementMatrix(elInfo,
452 453 454
				  elementMatrix, 
				  *factorIt ? **factorIt : 1.0);
      }
455
    }        
456

457
    addElementMatrix(factor, *elementMatrix, bound);   
458 459
  }

460 461 462 463 464 465 466 467
  void DOFMatrix::assemble(double factor, ElInfo *elInfo, const BoundaryType *bound,
			   Operator *op)
  {
      FUNCNAME("DOFMatrix::assemble()");

      TEST_EXIT_DBG(op)("No operator!\n");

      op->getAssembler(omp_get_thread_num())->initElementMatrix(elementMatrix, elInfo);
Thomas Witkowski's avatar
Thomas Witkowski committed
468
      op->getElementMatrix(elInfo, elementMatrix, factor);
469 470 471
      addElementMatrix(factor, *elementMatrix, bound);
  }

472 473 474
  void DOFMatrix::assemble(double factor, 
			   ElInfo *rowElInfo, ElInfo *colElInfo,
			   ElInfo *smallElInfo, ElInfo *largeElInfo,
475 476 477 478 479 480 481 482 483
			   const BoundaryType *bound, Operator *op)
  {
    FUNCNAME("DOFMatrix::assemble()");

    if (!op && operators.size() == 0) {
      return;
    }

    Operator *operat = op ? op : operators[0];
484 485
    operat->getAssembler(omp_get_thread_num())->
      initElementMatrix(elementMatrix, rowElInfo, colElInfo);
486 487
    
    if (op) {
Thomas Witkowski's avatar
Thomas Witkowski committed
488
      op->getElementMatrix(rowElInfo, colElInfo, smallElInfo, largeElInfo,
489
			   elementMatrix);
490
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
491 492 493
      std::vector<Operator*>::iterator it = operators.begin();
      std::vector<double*>::iterator factorIt = operatorFactor.begin();
      for (; it != operators.end(); ++it, ++factorIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531
	(*it)->getElementMatrix(rowElInfo, colElInfo,
				smallElInfo, largeElInfo,
				elementMatrix, 
				*factorIt ? **factorIt : 1.0);	
      }      
    }

    addElementMatrix(factor, *elementMatrix, bound);       
  }

  void DOFMatrix::assemble2(double factor, 
			    ElInfo *mainElInfo, ElInfo *auxElInfo,
			    ElInfo *smallElInfo, ElInfo *largeElInfo,			    
			   const BoundaryType *bound, Operator *op)
  {
    FUNCNAME("DOFMatrix::assemble2()");

    if (!op && operators.size() == 0) {
      return;
    }

    Operator *operat = op ? op : operators[0];
    operat->getAssembler(omp_get_thread_num())->
      initElementMatrix(elementMatrix, mainElInfo);
    
    if (op) {
      ERROR_EXIT("TODO");
//       op->getElementMatrix(rowElInfo, colElInfo, 
// 			   smallElInfo, largeElInfo,
// 			   elementMatrix);
    } else {
      std::vector<Operator*>::iterator it;
      std::vector<double*>::iterator factorIt;
      for (it = operators.begin(), factorIt = operatorFactor.begin();	
	   it != operators.end(); 
	   ++it, ++factorIt) {
	if ((*it)->getNeedDualTraverse()) {
	  (*it)->getElementMatrix(mainElInfo, auxElInfo,
532 533 534 535
				  smallElInfo, largeElInfo,
				  elementMatrix, 
				  *factorIt ? **factorIt : 1.0);
	}
536 537
      }      
    }
538

539 540 541
    addElementMatrix(factor, *elementMatrix, bound);       
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
542 543 544 545 546 547 548 549 550 551
  void DOFMatrix::finishAssembling()
  {
    // call the operatos cleanup procedures
    for (std::vector<Operator*>::iterator it = operators.begin();
	 it != operators.end();
	 ++it) {     
      (*it)->finishAssembling();
    }
  }

552 553 554
  Flag DOFMatrix::getAssembleFlag()
  {
    Flag fillFlag(0);
555
    std::vector<Operator*>::iterator op;
556 557 558
    for(op = operators.begin(); op != operators.end(); ++op) {
      fillFlag |= (*op)->getFillFlag();
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
559

560 561 562 563 564 565 566 567 568 569 570 571
    return fillFlag;
  }

  void DOFMatrix::mm(MatrixTranspose aTranspose,
		     DOFMatrix& a, 
		     MatrixTranspose bTranspose,
		     DOFMatrix& b)
  {
    FUNCNAME("DOFMatrix::mm()");

    WARNING("implementation not finished!!!\n");

572
    TEST_EXIT_DBG(a.getColFESpace() == b.getRowFESpace())
573
      ("a.colFESpace != b.rowFESpace\n");
574
    TEST_EXIT_DBG(rowFESpace == a.getRowFESpace())
575
      ("rowFESpace != a.rowFESpace\n");  
576
    TEST_EXIT_DBG(colFESpace == b.getColFESpace())
577 578 579 580
      ("colFESpace != b.colFESpace\n");  

    clear();

581
    int i, j;
582

583
    if (aTranspose == NoTranspose && bTranspose == NoTranspose) {
584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600
      int cols  = b.getColFESpace()->getAdmin()->getUsedSize();
      DOFMatrix::Iterator rowIterator(this, USED_DOFS);

      // for every row ...
      for(rowIterator.reset(); !rowIterator.end(); ++rowIterator) {
	int rowIndex = rowIterator.getDOFIndex();
	// and every column of result   
	for(i=0; i < cols; i++) {
	  double entry = 0.0;
	  // for every entry in a[i] ...
	  for(j=0; j < static_cast<int>( a[rowIndex].size()); j++) {
	    int logIndex = a[rowIndex][j].col;
	    int physIndex = b.logToPhysIndex(logIndex, i);
	    if(physIndex != -1) {
	      entry += a[rowIndex][j].entry * b[logIndex][physIndex].entry;
	    }
	  }
601
	  if (entry != 0.0) {
602 603 604 605 606 607 608 609 610 611 612
	    addSparseDOFEntry(1.0, rowIndex, i, entry);
	  }
	}
      }
    } else if(aTranspose == Transpose && bTranspose == NoTranspose) {
      DOFMatrix::Iterator aIterator(&a, USED_DOFS);
      DOFMatrix::Iterator bIterator(&b, USED_DOFS);
      for(aIterator.reset(), bIterator.reset(); 
	  !aIterator.end(); 
	  ++aIterator, ++bIterator) 
	{
613 614
	  std::vector<MatEntry>::const_iterator aRowIt;
	  std::vector<MatEntry>::const_iterator bRowIt;
615 616 617 618 619 620 621 622 623 624 625
	  for(aRowIt = aIterator->begin(); aRowIt != aIterator->end(); ++aRowIt) {
	    int aCol = aRowIt->col;
	    if(aCol == UNUSED_ENTRY) continue;
	    if(aCol == NO_MORE_ENTRIES) break;
	    for(bRowIt = bIterator->begin(); bRowIt !=bIterator->end(); ++bRowIt) {
	      int bCol = bRowIt->col;
	      if(bCol == UNUSED_ENTRY) continue;
	      if(bCol == NO_MORE_ENTRIES) break;

	      double entry = aRowIt->entry * bRowIt->entry;

626
	      if (entry != 0.0) {
627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642
		addSparseDOFEntry(1.0, aCol, bCol, entry);	    
	      }
	    }
	  }
	}
    } else if(aTranspose == NoTranspose && bTranspose == Transpose) {
      ERROR_EXIT("not yet\n");
    } else if(aTranspose == Transpose && bTranspose == Transpose) {
      ERROR_EXIT("not yet\n");
    }
  }

  void DOFMatrix::axpy(double a, 
		       const DOFMatrix& x,
		       const DOFMatrix& y)
  {
643 644 645 646
    FUNCNAME("DOFMatrix::axpy()");

    TEST_EXIT_DBG(x.getRowFESpace() == y.getRowFESpace() &&
		  rowFESpace == x.getRowFESpace())
647
      ("row fe-spaces not equal\n");
648 649
    TEST_EXIT_DBG(x.getColFESpace() == y.getColFESpace() &&
		  colFESpace == x.getColFESpace())
650
      ("col fe-spaces not equal\n");
651
    
652 653 654 655 656 657 658 659 660 661 662 663 664 665
    DOFMatrix::Iterator rowIterator(this, USED_DOFS);
    DOFMatrix::Iterator xIterator(const_cast<DOFMatrix*>(&x), USED_DOFS);
    DOFMatrix::Iterator yIterator(const_cast<DOFMatrix*>(&y), USED_DOFS);

    int i, rowIndex, colIndex;

    for(rowIterator.reset(), xIterator.reset(), yIterator.reset();
	!rowIterator.end();
	++rowIterator, ++xIterator, ++yIterator)
      {
	rowIndex = rowIterator.getDOFIndex();
	// add x contributions to this row
	for(i=0; i < static_cast<int>((*xIterator).size()); i++) {
	  colIndex = (*xIterator)[i].col;
666
	  if (colIndex >= 0) {
667 668 669 670 671 672
	    addSparseDOFEntry(a, rowIndex, colIndex, (*xIterator)[i].entry);
	  }
	}
	// add y contributions to this row
	for(i=0; i < static_cast<int>((*yIterator).size()); i++) {
	  colIndex = (*yIterator)[i].col;
673
	  if (colIndex >= 0) {
674 675 676 677 678 679 680 681 682
	    addSparseDOFEntry(1.0, rowIndex, colIndex, (*yIterator)[i].entry);
	  }
	}
      }
  }

  void DOFMatrix::scal(double b) 
  {
    DOFMatrix::Iterator rowIterator(this, USED_DOFS);
Thomas Witkowski's avatar
Thomas Witkowski committed
683 684 685
    for (rowIterator.reset(); !rowIterator.end(); ++rowIterator) {
      for (int i = 0; i < static_cast<int>((*rowIterator).size()); i++) {
	if ((*rowIterator)[i].col >= 0) {
686 687 688 689 690 691
	  (*rowIterator)[i].entry *= b;
	}
      }
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
692 693 694 695 696 697 698
  void DOFMatrix::addOperator(Operator *op, double* factor, double* estFactor) 
  { 
    operators.push_back(op);
    operatorFactor.push_back(factor);
    operatorEstFactor.push_back(estFactor);
  }

699 700 701 702 703
  void DOFMatrix::copy(const DOFMatrix& rhs) 
  {
    clear();
    DOFMatrix::Iterator rhsIterator(const_cast<DOFMatrix*>(&rhs), USED_DOFS);
    DOFMatrix::Iterator thisIterator(this, USED_DOFS);
704 705 706
    std::vector<MatEntry>::const_iterator colIt;
    std::vector<MatEntry>::const_iterator colBegin;
    std::vector<MatEntry>::const_iterator colEnd;
707 708 709 710 711 712 713 714 715 716 717 718 719 720 721
    for(rhsIterator.reset(), thisIterator.reset();
	!rhsIterator.end(); 
	++rhsIterator, ++thisIterator) 
      {
	colBegin = rhsIterator->begin();
	colEnd = rhsIterator->end();
	for(colIt = colBegin; colIt != colEnd; ++colIt) {
	  MatEntry matEntry;
	  matEntry.col = colIt->col;
	  matEntry.entry = colIt->entry;
	  thisIterator->push_back(matEntry);
	}
      }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
722 723 724 725 726 727 728 729 730 731 732 733 734
  void DOFMatrix::removeRowsWithDBC(std::set<int> *rows)
  {
    for (std::set<int>::iterator it = rows->begin();
	 it != rows->end(); 
	 ++it) {
      if (coupleMatrix) {
	matrix[*it].resize(0);
      } else {
	matrix[*it].resize(1);
	matrix[*it][0].col = *it;
	matrix[*it][0].entry = 1.0;
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
735 736

    rows->clear();
Thomas Witkowski's avatar
Thomas Witkowski committed
737 738
  }

739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755
  int DOFMatrix::getNumCols() const {
    int max = 0;
    int rowSize = static_cast<int>(matrix.size());
    for (int i = 0; i < rowSize; i++) {
      int colSize = static_cast<int>(matrix[i].size());
      for (int j = 0; j < colSize; j++) {
	if (matrix[i][j].col > max) {
	  max = matrix[i][j].col;
	}
      }
    }
    
    // Add one to the maximum value, because the indeces start with 0.
    return max + 1;
  }


756 757
  void DOFMatrix::createPictureFile(const char* filename, int dim)
  {
758 759
    png_structp png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, 
						  NULL, NULL, NULL);
760 761 762 763 764 765 766 767

    if (!png_ptr)
       return;

    png_bytep rowPointers[dim];
    for (int i = 0; i < dim; i++) {
      rowPointers[i] = (png_byte*)png_malloc(png_ptr, dim);

768
      for (int j = 0; j < dim; j++)
769 770 771 772 773 774 775 776
	rowPointers[i][j] = 255;
    }

    double scalFactor = static_cast<double>(dim) / static_cast<double>(matrix.size());

    for (int i = 0; i < static_cast<int>(matrix.size()); i++) {    
      int pi = static_cast<int>(static_cast<double>(i) * scalFactor);

777
      TEST_EXIT_DBG((pi >= 0) && (pi < dim))("PI");
778 779 780 781 782

      for (int j = 0; j < static_cast<int>(matrix[i].size()); j++) {
	
	int pj = static_cast<int>(static_cast<double>(matrix[i][j].col) * scalFactor);

783
	TEST_EXIT_DBG((pj >= 0) && (pj < dim))("PJ");
784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819

	rowPointers[pi][pj] = 0;
      }
    }
   
    FILE *fp = fopen(filename, "wb");
    TEST_EXIT(fp)("Cannot open file for writing matrix picture file!\n");

    png_infop info_ptr = png_create_info_struct(png_ptr);

    if (!info_ptr) {
       png_destroy_write_struct(&png_ptr, (png_infopp)NULL);
       return;
    }

    png_init_io(png_ptr, fp);

    png_set_IHDR(png_ptr, info_ptr, dim, dim, 8,
		 PNG_COLOR_TYPE_GRAY, 
		 PNG_INTERLACE_NONE,
		 PNG_COMPRESSION_TYPE_DEFAULT,
		 PNG_FILTER_TYPE_DEFAULT);

    png_set_rows(png_ptr, info_ptr, rowPointers);

    png_write_png(png_ptr, info_ptr, PNG_TRANSFORM_IDENTITY, NULL);

    png_destroy_write_struct(&png_ptr, &info_ptr);

    fclose(fp);
  }


  int DOFMatrix::memsize() 
  {   
    int sizeDOFMatrix = sizeof(DOFMatrix);   
820
    int sizeMatrix = sizeof(std::vector<MatrixRow>);
821 822 823 824 825 826 827 828
    for (int i = 0; i < static_cast<int>(matrix.size()); i++) {
      sizeMatrix += sizeof(MatrixRow) + matrix[i].size() * sizeof(MatEntry);
    }
    
    return sizeDOFMatrix + sizeMatrix;
  }


829
  double norm(std::vector<MatEntry> *row)
830 831
  {
    double result = 0.0;
832
    std::vector<MatEntry>::iterator it;
833 834 835 836 837 838 839
    for (it = row->begin(); it < row->end(); ++it) {
      result += (*it).entry * (*it).entry;
    }

    return(sqrt(result));
  }

840
  double min(std::vector<MatEntry> *row)
841 842 843 844 845 846
  {
    double result = 0.0;
    if (row->size() > 0) {
      result = (*row)[0].entry;
    }

847
    std::vector<MatEntry>::iterator it;
848 849 850 851 852 853 854 855 856
    for (it = row->begin(); it < row->end(); ++it) {
      if ((*it).entry < result) {
	result = (*it).entry;
      }
    }

    return(result);
  }

857
  double max(std::vector<MatEntry> *row)
858 859 860 861 862 863
  {
    double result = 0.0;
    if (row->size() > 0) {
      result = (*row)[0].entry;
    }

864
    std::vector<MatEntry>::iterator it;
865 866 867 868 869 870 871 872 873
    for (it = row->begin(); it < row->end(); ++it) {
      if ((*it).entry > result) {
	result = (*it).entry;
      }
    }

    return(result);
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916
  void addDOFMatrix(DOFMatrix *result, const DOFMatrix *a, const DOFMatrix *b)
  {
    result->clear();

    DOFMatrix::Iterator resultIterator(result, USED_DOFS);
    DOFMatrix::Iterator aIterator(const_cast<DOFMatrix*>(a), USED_DOFS);
    DOFMatrix::Iterator bIterator(const_cast<DOFMatrix*>(b), USED_DOFS);

    for (resultIterator.reset(), aIterator.reset(); 
	 !aIterator.end(); 
	 ++resultIterator, ++aIterator) {

      *resultIterator = *aIterator;
    }      

    for (resultIterator.reset(), bIterator.reset(); 
	 !bIterator.end(); 
	 ++resultIterator, ++bIterator) {
      std::vector<MatEntry>::iterator resultRowIt;
      std::vector<MatEntry>::const_iterator bRowIt;

      for (bRowIt = bIterator->begin();
	   bRowIt != bIterator->end();
	   ++bRowIt) {
	bool added = false;
	for (resultRowIt = resultIterator->begin();
	     resultRowIt != resultIterator->end();
	     ++resultRowIt) {
	  if (bRowIt->col == resultRowIt->col) {
	    resultRowIt->entry += bRowIt->entry;
	    
	    added = true;
	    break;
	  }
	}

	if (!added) {
	  resultIterator->push_back(*bRowIt);
	}
      }	         
    }      
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948
  void addDOFMatrix(DOFMatrix *result, const DOFMatrix *a)
  {
    DOFMatrix::Iterator resultIterator(result, USED_DOFS);
    DOFMatrix::Iterator aIterator(const_cast<DOFMatrix*>(a), USED_DOFS);

    for (resultIterator.reset(), aIterator.reset(); 
	 !aIterator.end(); 
	 ++resultIterator, ++aIterator) {
      std::vector<MatEntry>::iterator resultRowIt;
      std::vector<MatEntry>::const_iterator aRowIt;

      for (aRowIt = aIterator->begin();
	   aRowIt != aIterator->end();
	   ++aRowIt) {
	bool added = false;
	for (resultRowIt = resultIterator->begin();
	     resultRowIt != resultIterator->end();
	     ++resultRowIt) {
	  if (aRowIt->col == resultRowIt->col) {
	    resultRowIt->entry += aRowIt->entry;
	    
	    added = true;
	    break;
	  }
	}

	if (!added) {
	  resultIterator->push_back(*aRowIt);
	}
      }	         
    }      
  }
949
}