PetscSolverGlobalMatrix.cc 32.9 KB
Newer Older
Thomas Witkowski's avatar
Thomas Witkowski committed
1 2 3 4 5 6 7 8 9 10 11
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.

12 13
#include <mpi.h>
#include "DirichletBC.h"
Thomas Witkowski's avatar
Thomas Witkowski committed
14 15 16 17 18 19
#include "parallel/PetscSolverGlobalMatrix.h"
#include "parallel/StdMpi.h"
#include "parallel/MpiHelper.h"

namespace AMDiS {

20 21 22 23 24 25 26 27 28 29 30

  PetscSolverGlobalMatrix::PetscSolverGlobalMatrix(string name)
    : PetscSolver(name),
      zeroStartVector(false),
      printMatInfo(false)
  {
    Parameters::get("parallel->use zero start vector", zeroStartVector);
    Parameters::get("parallel->print matrix info", printMatInfo);
  }


Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
31
  void PetscSolverGlobalMatrix::fillPetscMatrix(Matrix<DOFMatrix*> *seqMat)
Thomas Witkowski's avatar
Thomas Witkowski committed
32 33 34
  {
    FUNCNAME("PetscSolverGlobalMatrix::fillPetscMatrix()");

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
35 36 37 38 39 40
    TEST_EXIT_DBG(meshDistributor)("No mesh distributor object defined!\n");
    TEST_EXIT_DBG(interiorMap)("No parallel mapping object defined!\n");
    TEST_EXIT_DBG(seqMat)("No DOF matrix defined!\n");
    
    double wtime = MPI::Wtime();

41
    createMatVec(*seqMat);
42

Thomas Witkowski's avatar
Thomas Witkowski committed
43
    if (coarseSpaceMap.size()) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
44
      fillPetscMatrixWithCoarseSpace(seqMat);
45 46 47
      return;
    }

48 49
    // === Create PETSc vector (solution and a temporary vector). ===

Thomas Witkowski's avatar
Thomas Witkowski committed
50 51 52 53 54 55
#if (DEBUG != 0)
    MSG("Fill petsc matrix 1 needed %.5f seconds\n", MPI::Wtime() - wtime);
#endif

    // === Transfer values from DOF matrices to the PETSc matrix. === 

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
56
    int nComponents = seqMat->getNumRows();
Thomas Witkowski's avatar
Thomas Witkowski committed
57 58
    for (int i = 0; i < nComponents; i++)
      for (int j = 0; j < nComponents; j++)
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
59 60
	if ((*seqMat)[i][j])
	  setDofMatrix((*seqMat)[i][j], i, j);
Thomas Witkowski's avatar
Thomas Witkowski committed
61 62 63 64 65

#if (DEBUG != 0)
    MSG("Fill petsc matrix 2 needed %.5f seconds\n", MPI::Wtime() - wtime);
#endif

66
    matAssembly();
67

Thomas Witkowski's avatar
Thomas Witkowski committed
68
    //    removeDirichletRows(seqMat);
69 70 71
    
    if (printMatInfo) {
      MatInfo matInfo;
72
      MatGetInfo(getMatInterior(), MAT_GLOBAL_SUM, &matInfo);
73 74 75 76 77 78 79
      MSG("Matrix info:\n");
      MSG("  memory usage: %e MB\n", matInfo.memory / (1024.0 * 1024.0));
      MSG("  mallocs: %d\n", static_cast<int>(matInfo.mallocs));
      MSG("  nz allocated: %d\n", static_cast<int>(matInfo.nz_allocated));
      MSG("  nz used: %d\n", static_cast<int>(matInfo.nz_used));
      MSG("  nz unneeded: %d\n", static_cast<int>(matInfo.nz_unneeded));
    }
80 81


82
    // === Init PETSc solver and preconditioner objects. ===
83

84 85
    initSolver(kspInterior);
    KSPGetPC(kspInterior, &pcInterior);
86
    initPreconditioner(pcInterior);
87

88

Thomas Witkowski's avatar
Thomas Witkowski committed
89 90 91
#if (DEBUG != 0)
    MSG("Fill petsc matrix 3 needed %.5f seconds\n", MPI::Wtime() - wtime);
#endif
92 93


94 95 96 97 98 99 100 101 102 103 104
    // === For debugging allow to write the matrix to a file. ===

    bool dbgWriteMatrix = false;
    Parameters::get("parallel->debug->write matrix", dbgWriteMatrix);
    if (dbgWriteMatrix) {
      PetscViewer matView;
      PetscViewerBinaryOpen(PETSC_COMM_WORLD, "mpi.mat",
			    FILE_MODE_WRITE, &matView);
      MatView(getMatInterior(), matView);
      PetscViewerDestroy(&matView);
    }
105 106 107
  }


Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
108
  void PetscSolverGlobalMatrix::fillPetscMatrixWithCoarseSpace(Matrix<DOFMatrix*> *seqMat)
109 110
  {
    FUNCNAME("PetscSolverGlobalMatrix::fillPetscMatrixWithCoarseSpace()");
111 112

    TEST_EXIT_DBG(interiorMap)("Should not happen!\n");
113 114
    TEST_EXIT_DBG(coarseSpaceMap.size() == seqMat->getSize())
      ("Wrong sizes %d %d\n", coarseSpaceMap.size(), seqMat->getSize());
115

116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134
    // === Prepare traverse of sequentially created matrices. ===

    using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
    namespace traits = mtl::traits;
    typedef DOFMatrix::base_matrix_type Matrix;

    typedef traits::range_generator<row, Matrix>::type cursor_type;
    typedef traits::range_generator<nz, cursor_type>::type icursor_type;

    vector<int> cols, colsOther;
    vector<double> values, valuesOther;
    cols.reserve(300);
    colsOther.reserve(300);
    values.reserve(300);
    valuesOther.reserve(300);

    // === Traverse all sequentially created matrices and add the values to ===
    // === the global PETSc matrices.                                       ===

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
135
    int nComponents = seqMat->getSize();
Thomas Witkowski's avatar
Thomas Witkowski committed
136 137 138 139 140
    for (int rowComponent = 0; rowComponent < nComponents; rowComponent++) {
      for (int colComponent = 0; colComponent < nComponents; colComponent++) {
	DOFMatrix* dofMat = (*seqMat)[rowComponent][colComponent];

	if (!dofMat)
141 142
	  continue;

Thomas Witkowski's avatar
Thomas Witkowski committed
143 144
	ParallelDofMapping *rowCoarseSpace = coarseSpaceMap[rowComponent];
	ParallelDofMapping *colCoarseSpace = coarseSpaceMap[colComponent];
Thomas Witkowski's avatar
Thomas Witkowski committed
145

Thomas Witkowski's avatar
Thomas Witkowski committed
146 147
	std::set<DegreeOfFreedom> &dirichletRows = dofMat->getDirichletRows();

Thomas Witkowski's avatar
Thomas Witkowski committed
148 149
	traits::col<Matrix>::type col(dofMat->getBaseMatrix());
	traits::const_value<Matrix>::type value(dofMat->getBaseMatrix());
150 151
	
	// Traverse all rows.
Thomas Witkowski's avatar
Thomas Witkowski committed
152 153
	for (cursor_type cursor = begin<row>(dofMat->getBaseMatrix()), 
	       cend = end<row>(dofMat->getBaseMatrix()); cursor != cend; ++cursor) {
154

Thomas Witkowski's avatar
Thomas Witkowski committed
155
	  bool isRowCoarse = isCoarseSpace(rowComponent, *cursor);
Thomas Witkowski's avatar
Thomas Witkowski committed
156 157 158 159 160 161 162 163

	  // For the case, this is a dirichlet row we have to check whether the 
	  // rank is also owner of this row DOF.
	  if (dirichletRows.count(*cursor)) {
	    if ((!isRowCoarse && !(*interiorMap)[rowComponent].isRankDof(*cursor)) ||
		(isRowCoarse && !(*rowCoarseSpace)[rowComponent].isRankDof(*cursor)))
	      continue;
	  }
164 165 166 167 168 169 170 171 172 173
  
	  cols.clear();
	  colsOther.clear();
	  values.clear();	  
	  valuesOther.clear();

	  // Traverse all columns.
	  for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); 
	       icursor != icend; ++icursor) {

Thomas Witkowski's avatar
Thomas Witkowski committed
174
	    bool isColCoarse = isCoarseSpace(colComponent, col(*icursor));
175

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
176
	    if (isColCoarse == false)
177
	      if ((*interiorMap)[colComponent].isSet(col(*icursor)) == false)
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
178 179
		continue;

Thomas Witkowski's avatar
Thomas Witkowski committed
180 181 182
	    if (isColCoarse == isRowCoarse) {
	      cols.push_back(col(*icursor));
	      values.push_back(value(*icursor));
183
	    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
184 185
	      colsOther.push_back(col(*icursor));
	      valuesOther.push_back(value(*icursor));
186 187 188 189 190 191
	    }
	  }  // for each nnz in row


	  // === Set matrix values. ===

Thomas Witkowski's avatar
Thomas Witkowski committed
192
	  if (isRowCoarse) {
Thomas Witkowski's avatar
Thomas Witkowski committed
193 194
	    for (unsigned int i = 0; i < cols.size(); i++)
	      cols[i] = colCoarseSpace->getMatIndex(colComponent, cols[i]);
195

Thomas Witkowski's avatar
Thomas Witkowski committed
196
	    int rowIndex = rowCoarseSpace->getMatIndex(rowComponent, *cursor);
Thomas Witkowski's avatar
Thomas Witkowski committed
197
	    MatSetValues(getMatCoarseByComponent(rowComponent, colComponent),
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
198
			 1, &rowIndex, cols.size(),
199 200 201
			 &(cols[0]), &(values[0]), ADD_VALUES);

	    if (colsOther.size()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
202 203 204
	      for (unsigned int i = 0; i < colsOther.size(); i++)
		colsOther[i] = 
		  interiorMap->getMatIndex(colComponent, colsOther[i]) + 
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
205
		  rStartInterior;	      
206

Thomas Witkowski's avatar
Thomas Witkowski committed
207
	      MatSetValues(getMatCoarseInteriorByComponent(rowComponent), 
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
208
			   1, &rowIndex, colsOther.size(),
209 210 211
 			   &(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
	    }
	  } else {
212
	    if ((*interiorMap)[rowComponent].isSet(*cursor) == false)
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
213 214
	      continue;

215
	    int localRowIndex = 
Thomas Witkowski's avatar
Thomas Witkowski committed
216 217 218
	      (subdomainLevel == 0 ? 
	       interiorMap->getLocalMatIndex(rowComponent, *cursor) :
	       interiorMap->getMatIndex(rowComponent, *cursor));
219

Thomas Witkowski's avatar
Thomas Witkowski committed
220
	    for (unsigned int i = 0; i < cols.size(); i++) {
221
	      if (subdomainLevel == 0)
Thomas Witkowski's avatar
Thomas Witkowski committed
222
		cols[i] = interiorMap->getLocalMatIndex(colComponent, cols[i]);
223
	      else
Thomas Witkowski's avatar
Thomas Witkowski committed
224
		cols[i] = interiorMap->getMatIndex(colComponent, cols[i]);
225 226
	    }
	    
227
  	    MatSetValues(getMatInterior(), 1, &localRowIndex, cols.size(),
228
  			 &(cols[0]), &(values[0]), ADD_VALUES);
Thomas Witkowski's avatar
Thomas Witkowski committed
229
	    
230
	    if (colsOther.size()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
231 232 233
	      int globalRowIndex = 
		interiorMap->getMatIndex(rowComponent, *cursor) + rStartInterior;
      
Thomas Witkowski's avatar
Thomas Witkowski committed
234 235 236
	      for (unsigned int i = 0; i < colsOther.size(); i++)
		colsOther[i] = 
		  colCoarseSpace->getMatIndex(colComponent, colsOther[i]);
237

Thomas Witkowski's avatar
Thomas Witkowski committed
238
  	      MatSetValues(getMatInteriorCoarseByComponent(colComponent), 
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
239
			   1, &globalRowIndex, colsOther.size(),
240 241 242 243 244 245 246
  			   &(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
	    }
	  }
	} 
      }
    }

247
    matAssembly();
248

Thomas Witkowski's avatar
Thomas Witkowski committed
249
    //    removeDirichletRows(seqMat);
250

251 252 253
    // === Create solver for the non primal (thus local) variables. ===

    KSPCreate(mpiCommLocal, &kspInterior);
254
    KSPSetOperators(kspInterior, getMatInterior(), getMatInterior(),
255
		    SAME_NONZERO_PATTERN);
256 257 258
    KSPSetOptionsPrefix(kspInterior, "interior_");
    KSPSetType(kspInterior, KSPPREONLY);
    KSPGetPC(kspInterior, &pcInterior);
259 260
    if (isSymmetric) {
      PCSetType(pcInterior, PCCHOLESKY);
261
      PCFactorSetMatSolverPackage(pcInterior, MATSOLVERMUMPS);
262 263 264 265 266 267 268
    } else {
      PCSetType(pcInterior, PCLU);
      if (subdomainLevel == 0)
	PCFactorSetMatSolverPackage(pcInterior, MATSOLVERUMFPACK);
      else
	PCFactorSetMatSolverPackage(pcInterior, MATSOLVERMUMPS);
    }
269
    KSPSetFromOptions(kspInterior);
270 271 272
  }


273 274 275
  void PetscSolverGlobalMatrix::fillPetscRhs(SystemVector *vec)
  {
    FUNCNAME("PetscSolverGlobalMatrix::fillPetscRhs()");
Thomas Witkowski's avatar
Thomas Witkowski committed
276

277
    TEST_EXIT_DBG(vec)("No DOF vector defined!\n");
278
    TEST_EXIT_DBG(interiorMap)("No parallel DOF map defined!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
279
    
Thomas Witkowski's avatar
Thomas Witkowski committed
280
    // === Transfer values from DOF vector to the PETSc vector. === 
Thomas Witkowski's avatar
Thomas Witkowski committed
281
    if (coarseSpaceMap.size()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
282
      for (int i = 0; i < vec->getSize(); i++)
283
	setDofVector(getVecRhsInterior(), 
Thomas Witkowski's avatar
Thomas Witkowski committed
284
		     getVecRhsCoarseByComponent(i), vec->getDOFVector(i), i);
Thomas Witkowski's avatar
Thomas Witkowski committed
285 286
    } else {
      for (int i = 0; i < vec->getSize(); i++)
287
	setDofVector(getVecRhsInterior(), vec->getDOFVector(i), i);
Thomas Witkowski's avatar
Thomas Witkowski committed
288 289
    }

290
    vecRhsAssembly();
291

Thomas Witkowski's avatar
Thomas Witkowski committed
292
    //    removeDirichletRows(vec);
293 294 295 296 297 298 299 300 301 302 303 304

    // === For debugging allow to write the rhs vector to a file. ===

    bool dbgWriteRhs = false;
    Parameters::get("parallel->debug->write rhs", dbgWriteRhs);
    if (dbgWriteRhs) {
      PetscViewer vecView;
      PetscViewerBinaryOpen(PETSC_COMM_WORLD, "mpi.vec",
			    FILE_MODE_WRITE, &vecView);
      VecView(getVecRhsInterior(), vecView);
      PetscViewerDestroy(&vecView);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
305 306 307
  }


308 309
  void PetscSolverGlobalMatrix::solvePetscMatrix(SystemVector &vec, 
						 AdaptInfo *adaptInfo)
Thomas Witkowski's avatar
Thomas Witkowski committed
310 311 312 313 314 315 316
  {
    FUNCNAME("PetscSolverGlobalMatrix::solvePetscMatrix()");

    int nComponents = vec.getSize();

    // === Set old solution to be initiual guess for PETSc solver. ===
    if (!zeroStartVector) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
317 318
      TEST_EXIT(coarseSpaceMap.size() == 0)("Not yet supported!\n");

319
      VecSet(getVecSolInterior(), 0.0);
Thomas Witkowski's avatar
Thomas Witkowski committed
320 321
      
      for (int i = 0; i < nComponents; i++)
322
	setDofVector(getVecSolInterior(), vec.getDOFVector(i), i, true);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
323

324
      vecSolAssembly();
Thomas Witkowski's avatar
Thomas Witkowski committed
325 326
    }

327

Thomas Witkowski's avatar
Thomas Witkowski committed
328 329
    MatNullSpace matNullspace;
    Vec nullspaceBasis;
330 331 332
    if (nullspace.size() > 0 || 
	hasConstantNullspace ||
	constNullspaceComponent.size() > 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
333
      TEST_EXIT_DBG(nullspace.size() <= 1)("Not yet implemented!\n");
334

335 336 337 338 339 340 341 342 343 344
      if (constNullspaceComponent.size() > 0) {
	nullspace.clear();
	SystemVector *basisVec = new SystemVector(vec);
	basisVec->set(0.0);
	for (unsigned int i = 0; i < constNullspaceComponent.size(); i++)
	  basisVec->getDOFVector(constNullspaceComponent[i])->set(1.0);

	nullspace.push_back(basisVec);
      } 

Thomas Witkowski's avatar
Thomas Witkowski committed
345
      if (nullspace.size() > 0) {
346
	VecDuplicate(getVecSolInterior(), &nullspaceBasis);
347
	setDofVector(nullspaceBasis, *(nullspace[0]), true);
Thomas Witkowski's avatar
Thomas Witkowski committed
348 349 350
	
	VecAssemblyBegin(nullspaceBasis);
	VecAssemblyEnd(nullspaceBasis);
351 352

	VecNormalize(nullspaceBasis, PETSC_NULL);
Thomas Witkowski's avatar
Thomas Witkowski committed
353 354 355 356
	
	MatNullSpaceCreate(mpiCommGlobal, (hasConstantNullspace ? PETSC_TRUE : PETSC_FALSE), 
			   1, &nullspaceBasis, &matNullspace);

357
	MatMult(getMatInterior(), nullspaceBasis, getVecSolInterior());
Thomas Witkowski's avatar
Thomas Witkowski committed
358
	PetscReal n;
359
	VecNorm(getVecSolInterior(), NORM_2, &n);
Thomas Witkowski's avatar
Thomas Witkowski committed
360 361 362 363
	MSG("NORM IS: %e\n", n);
      } else {
	MatNullSpaceCreate(mpiCommGlobal, PETSC_TRUE, 0, PETSC_NULL, &matNullspace);
      }
364

Thomas Witkowski's avatar
Thomas Witkowski committed
365 366 367
      MSG("NULLSPACE IS NOT REMOVED!\n");
      // MatSetNullSpace(getMatInterior(), matNullspace);
      // KSPSetNullSpace(kspInterior, matNullspace);
368

Thomas Witkowski's avatar
Thomas Witkowski committed
369
      // === Remove null space, if requested. ===
370

Thomas Witkowski's avatar
Thomas Witkowski committed
371
      if (removeRhsNullspace) {
Thomas Witkowski's avatar
Thomas Witkowski committed
372
	TEST_EXIT_DBG(coarseSpaceMap.empty())("Not supported!\n");
373 374

	MSG("Remove nullspace from rhs vector.\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
375
	
376
	MatNullSpaceRemove(matNullspace, getVecRhsInterior(), PETSC_NULL);
Thomas Witkowski's avatar
Thomas Witkowski committed
377 378 379 380
      }
    } else {
      TEST_EXIT(removeRhsNullspace == false)
	("No nullspace provided that should be removed from rhs!\n");
381 382
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
383
    // PETSc.
384
    solve(getVecRhsInterior(), getVecSolInterior());
Thomas Witkowski's avatar
Thomas Witkowski committed
385

386 387

    if (nullspace.size() > 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
388
      MatNullSpaceDestroy(&matNullspace);
389 390 391 392
      VecDestroy(&nullspaceBasis);
    }


Thomas Witkowski's avatar
Thomas Witkowski committed
393 394
    // === Transfere values from PETSc's solution vectors to the DOF vectors. ===
    PetscScalar *vecPointer;
395
    VecGetArray(getVecSolInterior(), &vecPointer);    
Thomas Witkowski's avatar
Thomas Witkowski committed
396

397
    int c = 0;
398 399
    for (int component = 0; component < nComponents; component++) {
      DOFVector<double> &dv = *(vec.getDOFVector(component));
400

401
      DofMap& d = (*interiorMap)[component].getMap();
402 403 404
      for (DofMap::iterator it = d.begin(); it != d.end(); ++it)
	if (it->second.local != -1)
	  dv[it->first] = vecPointer[c++];
Thomas Witkowski's avatar
Thomas Witkowski committed
405 406
    }

407
    VecRestoreArray(getVecSolInterior(), &vecPointer);
Thomas Witkowski's avatar
Thomas Witkowski committed
408 409 410 411

    // === Synchronize DOFs at common DOFs, i.e., DOFs that correspond to ===  
    // === more than one partition.                                       ===  
    meshDistributor->synchVector(vec);
412 413 414
  }


415 416 417 418
  void PetscSolverGlobalMatrix::solveGlobal(Vec &rhs, Vec &sol)
  {
    FUNCNAME("PetscSolverGlobalMatrix::solveGlobal()");

419 420 421
    double wtime = MPI::Wtime();
    double t0 = 0.0, t1 = 0.0;

422 423
    Vec tmp;
    if (mpiCommLocal.Get_size() == 1)
Thomas Witkowski's avatar
Thomas Witkowski committed
424
      interiorMap->createLocalVec(tmp);
425
    else
Thomas Witkowski's avatar
Thomas Witkowski committed
426
      interiorMap->createVec(tmp);
427 428 429 430 431 432 433 434 435 436 437

    PetscScalar *tmpValues, *rhsValues;
    VecGetArray(tmp, &tmpValues);
    VecGetArray(rhs, &rhsValues);

    for (int i = 0; i < interiorMap->getRankDofs(); i++)
      tmpValues[i] = rhsValues[i];

    VecRestoreArray(rhs, &rhsValues);
    VecRestoreArray(tmp, &tmpValues);

438 439 440
    t0 = MPI::Wtime() - wtime;

    wtime = MPI::Wtime();
441
    KSPSolve(kspInterior, tmp, tmp);
442 443 444
    t1 = MPI::Wtime() - wtime;

    wtime = MPI::Wtime();
445 446 447 448 449 450 451 452 453 454 455

    VecGetArray(tmp, &tmpValues);
    VecGetArray(sol, &rhsValues);

    for (int i = 0; i < interiorMap->getRankDofs(); i++) 
      rhsValues[i] = tmpValues[i];

    VecRestoreArray(sol, &rhsValues);
    VecRestoreArray(tmp, &tmpValues);

    VecDestroy(&tmp);
456 457
    t0 += MPI::Wtime() - wtime;

458
    //    MSG("TIMEING: %.5f %.5f\n", t0, t1);
459 460 461
  }


462 463 464 465
  void PetscSolverGlobalMatrix::destroyMatrixData()
  {
    FUNCNAME("PetscSolverGlobalMatrix::destroyMatrixData()");

466
    matDestroy();
Thomas Witkowski's avatar
Thomas Witkowski committed
467

468 469 470
    exitPreconditioner(pcInterior);

    exitSolver(kspInterior);
Thomas Witkowski's avatar
Thomas Witkowski committed
471 472 473
  }


474 475 476 477
  void PetscSolverGlobalMatrix::destroyVectorData()
  {
    FUNCNAME("PetscSolverGlobalMatrix::destroyVectorData()");

478
    vecDestroy();
479 480 481
  }


482
  void PetscSolverGlobalMatrix::removeDirichletRows(Matrix<DOFMatrix*> *seqMat)					       
483 484 485 486 487 488
  {
    FUNCNAME("PetscSolverGlobalMatrix::removeDirichletRows()");

    if (!handleDirichletRows)
      return;

489 490 491 492
    vector<int> dRowsInterior, dRowsCoarse;
    vector<int> dColsInterior, dColsCoarse;
    vector<double> dValuesInterior, dValuesCoarse;

493 494
    int nComponents = seqMat->getSize();

495
    for (int rowComp = 0; rowComp < nComponents; rowComp++) {
496 497 498
      bool dirichletRow = false;
      int dirichletMainCol = -1;

499 500
      for (int colComp = 0; colComp < nComponents; colComp++) {
	DOFMatrix *dofMat = (*seqMat)[rowComp][colComp];
501 502 503 504 505 506 507 508 509 510
	if (!dofMat)
	  continue;

	BoundaryIndexMap& bounds = 
	  const_cast<BoundaryIndexMap&>(dofMat->getBoundaryManager()->getBoundaryConditionMap());

	for (BoundaryIndexMap::iterator bIt = bounds.begin(); bIt != bounds.end(); ++bIt) {
	  if (bIt->second && bIt->second->isDirichlet()) {
	    dirichletRow = true;
	    if ((dynamic_cast<DirichletBC*>(bIt->second))->applyBoundaryCondition()) {
511
	      dirichletMainCol = colComp;
512
	      break;
513 514 515 516 517 518
	    } 
	  }
	}
      }

      if (!dirichletRow)
519 520
	continue;

521
      DOFMatrix *dofMat = (*seqMat)[rowComp][dirichletMainCol];
522 523
      TEST_EXIT(dofMat->getRowFeSpace() == dofMat->getColFeSpace())
	("I have to think about this scenario! Really possible?\n");
524

525

526 527 528
      std::set<DegreeOfFreedom> &dRows = dofMat->getDirichletRows();
      for (std::set<DegreeOfFreedom>::iterator dofIt = dRows.begin();
	   dofIt != dRows.end(); ++dofIt) {
529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569
	if (hasCoarseSpace()) {
	  bool isRowCoarse = isCoarseSpace(rowComp, *dofIt);
	  bool isColCoarse = isCoarseSpace(dirichletMainCol, *dofIt);
	  TEST_EXIT(isRowCoarse == isColCoarse)
	    ("Really possible? So reimplement AMDiS from the scratch!\n");

	  if (isRowCoarse) {
	    ParallelDofMapping *rowCoarseSpace = coarseSpaceMap[rowComp];
	    ParallelDofMapping *colCoarseSpace = coarseSpaceMap[dirichletMainCol];

	    MultiIndex multiIndex;
	    if ((*rowCoarseSpace)[rowComp].find(*dofIt, multiIndex) &&
		(*rowCoarseSpace)[rowComp].isRankDof(*dofIt)) {
	      int rowIndex = rowCoarseSpace->getLocalMatIndex(rowComp, *dofIt);
	      int colIndex = colCoarseSpace->getLocalMatIndex(dirichletMainCol, *dofIt);
	      dRowsCoarse.push_back(rowIndex);
	      dColsCoarse.push_back(colIndex);
	      dValuesCoarse.push_back(1.0);
	    }
	  } else {

	    MultiIndex multiIndex;
	    if ((*interiorMap)[rowComp].find(*dofIt, multiIndex) &&
		(*interiorMap)[rowComp].isRankDof(*dofIt)) {
	      int rowIndex = interiorMap->getLocalMatIndex(rowComp, *dofIt);
	      int colIndex = interiorMap->getLocalMatIndex(dirichletMainCol, *dofIt);
	      dRowsInterior.push_back(rowIndex);
	      dColsInterior.push_back(colIndex);
	      dValuesInterior.push_back(1.0);
	    }
	  }
	} else {
	  MultiIndex multiIndex;
	  if ((*interiorMap)[rowComp].find(*dofIt, multiIndex) &&
	      (*interiorMap)[rowComp].isRankDof(*dofIt)) {
	    int rowIndex = interiorMap->getMatIndex(rowComp, multiIndex.global);
	    int colIndex = interiorMap->getMatIndex(dirichletMainCol, multiIndex.global);
	    dRowsInterior.push_back(rowIndex);
	    dColsInterior.push_back(colIndex);
	    dValuesInterior.push_back(1.0);
	  }
570
	}
571 572 573
      }
    }

574 575 576 577 578 579 580 581 582 583 584 585
    {
      Mat mpiMat = getMatInterior();
      MatZeroRows(mpiMat, dRowsInterior.size(), &(dRowsInterior[0]), 0.0,
		  PETSC_NULL, PETSC_NULL);
      
      for (int i = 0; i < static_cast<int>(dRowsInterior.size()); i++)
	MatSetValue(mpiMat, dRowsInterior[i], dColsInterior[i], 
		    dValuesInterior[i], INSERT_VALUES);    
      
      MatAssemblyBegin(mpiMat, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(mpiMat, MAT_FINAL_ASSEMBLY);
    }
586

587
    if (hasCoarseSpace()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
588 589 590 591 592 593 594 595 596
      MatZeroRows(getMatInteriorCoarse(0), 
		  dRowsInterior.size(), &(dRowsInterior[0]), 0.0, PETSC_NULL, PETSC_NULL);
      MatZeroRows(getMatInteriorCoarse(1), 
		  dRowsInterior.size(), &(dRowsInterior[0]), 0.0, PETSC_NULL, PETSC_NULL);


      MatZeroRows(getMatCoarseInterior(0), 
		  dRowsCoarse.size(), &(dRowsCoarse[0]), 0.0, PETSC_NULL, PETSC_NULL);

597 598 599 600 601 602 603 604 605 606 607
      Mat mpiMat = getMatCoarse();
      MatZeroRows(mpiMat, dRowsCoarse.size(), &(dRowsCoarse[0]), 0.0,
		  PETSC_NULL, PETSC_NULL);
      
      for (int i = 0; i < static_cast<int>(dRowsCoarse.size()); i++)
	MatSetValue(mpiMat, dRowsCoarse[i], dColsCoarse[i], 
		    dValuesCoarse[i], INSERT_VALUES);    
      
      MatAssemblyBegin(mpiMat, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(mpiMat, MAT_FINAL_ASSEMBLY);
    }
608 609 610
  }

  
611
  void PetscSolverGlobalMatrix::removeDirichletRows(SystemVector *seqVec)
612 613 614 615 616 617
  {
    FUNCNAME("PetscSolverGlobalMatrix::removeDirichletRows()");

    if (!handleDirichletRows)
      return;

618 619
    int cInterior = 0;
    int cCoarse = 0;
620

621 622 623
    int nComponents = seqVec->getSize();
    for (int component = 0; component < nComponents; component++) {
      DOFVector<double> *dofVec = seqVec->getDOFVector(component);
624
      map<DegreeOfFreedom, double>& dValues = dofVec->getDirichletValues();
625

626 627
      for (map<DegreeOfFreedom, double>::iterator dofIt = dValues.begin();
	   dofIt != dValues.end(); ++dofIt) {
628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661
	if (hasCoarseSpace()) {
	  if (isCoarseSpace(component, dofIt->first)) {
	    ParallelDofMapping *rowCoarseSpace = coarseSpaceMap[component];
	    MultiIndex multiIndex;
	    if ((*rowCoarseSpace)[component].find(dofIt->first, multiIndex) &&
		(*rowCoarseSpace)[component].isRankDof(dofIt->first)) {
	      if (dofIt->second > 0.5)
		cCoarse++;
	      VecSetValue(getVecRhsCoarse(), 
			  rowCoarseSpace->getMatIndex(component, multiIndex.global),
			  dofIt->second, INSERT_VALUES);
	    }
	  } else {
	    MultiIndex multiIndex;
	    if ((*interiorMap)[component].find(dofIt->first, multiIndex) &&
		(*interiorMap)[component].isRankDof(dofIt->first)) {
	      if (dofIt->second > 0.5)
		cInterior++;
	      VecSetValue(getVecRhsInterior(), 
			  interiorMap->getLocalMatIndex(component, dofIt->first),
			  dofIt->second, INSERT_VALUES);
	    }
	  }
	} else {
	  MultiIndex multiIndex;
	  if ((*interiorMap)[component].find(dofIt->first, multiIndex) &&
	      (*interiorMap)[component].isRankDof(dofIt->first)) {
	    if (dofIt->second > 0.5)
	      cInterior++;
	    VecSetValue(getVecRhsInterior(), 
			interiorMap->getMatIndex(component, multiIndex.global),
			dofIt->second, INSERT_VALUES);
	  }
	}
662 663 664
      }
    }

665 666 667 668 669 670 671 672 673 674 675
    VecAssemblyBegin(getVecRhsInterior());
    VecAssemblyEnd(getVecRhsInterior());

    if (hasCoarseSpace()) {
      VecAssemblyBegin(getVecRhsCoarse());
      VecAssemblyEnd(getVecRhsCoarse());
    }

    mpi::globalAdd(cInterior);
    mpi::globalAdd(cCoarse);
    MSG("WORKED ON DIRICHLET DOFS: %d %d\n", cInterior, cCoarse);
676 677 678
  }


679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697
  void PetscSolverGlobalMatrix::createFieldSplit(PC pc)
  {
    FUNCNAME("PetscSolverGlobalMatrix::createFieldSplit()");

    vector<string> isNames;
    Parameters::get("parallel->solver->is blocks", isNames);

    int nBlocks = isNames.size();
    if (nBlocks == 0)
      return;

    for (int i = 0; i < nBlocks; i++) {
      MSG("Create for block %s\n", isNames[i].c_str());

      vector<int> blockComponents;
      Parameters::get("parallel->solver->is block " + lexical_cast<string>(i),
		      blockComponents);
      int nComponents = static_cast<int>(blockComponents.size());

698
      TEST_EXIT(nComponents > 0)("No IS block for block %d defined!\n", i);
699 700 701 702 703 704 705 706

      // Check if blocks are continous
      for (int j = 0; j < nComponents; j++) {
	TEST_EXIT(blockComponents[j] == blockComponents[0] + j)
	  ("Does not yet support not continous IS blocks! Block %s\n", 
	   isNames[i].c_str());
      }

707
      createFieldSplit(pc, isNames[i], blockComponents);
708 709 710 711
    }
  }


712 713 714 715 716 717 718 719
  void PetscSolverGlobalMatrix::createFieldSplit(PC pc, 
						 string splitName, 
						 vector<int> &components)
  {
    FUNCNAME("PetscSolverGlobalMatrix::createFieldSplit()");

    IS is;
    interiorMap->createIndexSet(is, components[0], components.size());
Thomas Witkowski's avatar
Thomas Witkowski committed
720
    PCFieldSplitSetIS(pc, splitName.c_str(), is);
721 722 723 724
    ISDestroy(&is);
  }


725
  void PetscSolverGlobalMatrix::initSolver(KSP &ksp)
726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750
  {
    FUNCNAME("PetscSolverGlobalMatrix::initSolver()");

    KSPCreate(mpiCommGlobal, &ksp);
    KSPSetOperators(ksp, getMatInterior(), getMatInterior(), 
		    SAME_NONZERO_PATTERN); 
    KSPSetTolerances(ksp, 0.0, 1e-8, PETSC_DEFAULT, PETSC_DEFAULT);
    KSPSetType(ksp, KSPBCGS);
    KSPSetOptionsPrefix(ksp, kspPrefix.c_str());
    KSPSetFromOptions(ksp);

    // Do not delete the solution vector, use it for the initial guess.
    if (!zeroStartVector)
      KSPSetInitialGuessNonzero(ksp, PETSC_TRUE);
  }  


  void PetscSolverGlobalMatrix::exitSolver(KSP ksp)
  {
    FUNCNAME("PetscSolverGlobalMatrix::exitSolver()");

    KSPDestroy(&ksp);
  }  


751 752 753 754 755 756 757 758 759 760 761 762 763 764 765
  void PetscSolverGlobalMatrix::initPreconditioner(PC pc)
  {
    FUNCNAME("PetscSolverGlobalMatrix::initPreconditioner()");

    PCSetFromOptions(pc);
    createFieldSplit(pc);
  }

 
  void PetscSolverGlobalMatrix::exitPreconditioner(PC pc)
  {
    FUNCNAME("PetscSolverGlobalMatrix::exitPreconditioner()");
  }


Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
766
  void PetscSolverGlobalMatrix::setDofMatrix(DOFMatrix* seqMat,
Thomas Witkowski's avatar
Thomas Witkowski committed
767
					     int rowComp, int colComp)
Thomas Witkowski's avatar
Thomas Witkowski committed
768 769 770
  {
    FUNCNAME("PetscSolverGlobalMatrix::setDofMatrix()");

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
771
    TEST_EXIT(seqMat)("No DOFMatrix!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
772 773

    using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
774
    namespace traits = mtl::traits;
Thomas Witkowski's avatar
Thomas Witkowski committed
775 776
    typedef DOFMatrix::base_matrix_type Matrix;

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
777 778
    traits::col<Matrix>::type col(seqMat->getBaseMatrix());
    traits::const_value<Matrix>::type value(seqMat->getBaseMatrix());
Thomas Witkowski's avatar
Thomas Witkowski committed
779 780 781 782 783 784 785 786 787 788 789

    typedef traits::range_generator<row, Matrix>::type cursor_type;
    typedef traits::range_generator<nz, cursor_type>::type icursor_type;

    vector<int> cols;
    vector<double> values;
    cols.reserve(300);
    values.reserve(300);
    
    vector<int> globalCols;

790 791
    // Get periodic mapping object
    PeriodicMap &perMap = meshDistributor->getPeriodicMap();
Thomas Witkowski's avatar
Thomas Witkowski committed
792 793
    std::set<DegreeOfFreedom> &dirichletRows = seqMat->getDirichletRows();
      
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
794 795
    const FiniteElemSpace *rowFe = seqMat->getRowFeSpace();
    const FiniteElemSpace *colFe = seqMat->getColFeSpace();
Thomas Witkowski's avatar
Thomas Witkowski committed
796 797

    // === Traverse all rows of the DOF matrix and insert row wise the values ===
Thomas Witkowski's avatar
Thomas Witkowski committed
798 799
    // === to the PETSc matrix.                                               ===

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
800 801
    for (cursor_type cursor = begin<row>(seqMat->getBaseMatrix()), 
	   cend = end<row>(seqMat->getBaseMatrix()); cursor != cend; ++cursor) {
Thomas Witkowski's avatar
Thomas Witkowski committed
802
      // Global index of the current row DOF.
803
      MultiIndex rowMultiIndex;
Thomas Witkowski's avatar
Thomas Witkowski committed
804
      if ((*interiorMap)[rowComp].find(*cursor, rowMultiIndex) == false)
805 806 807
	continue;

      int globalRowDof = rowMultiIndex.global;
808

Thomas Witkowski's avatar
Thomas Witkowski committed
809
      // Test if the current row DOF is a periodic DOF.
810
      bool periodicRow = perMap.isPeriodic(rowFe, globalRowDof);
811

Thomas Witkowski's avatar
Thomas Witkowski committed
812 813 814 815
      // Dirichlet rows can be set only be the owner ranks.
      if (dirichletRows.count(*cursor) && !((*interiorMap)[rowComp].isRankDof(*cursor)))
	continue;

Thomas Witkowski's avatar
Thomas Witkowski committed
816 817 818
      if (!periodicRow) {
	// === Row DOF index is not periodic. ===

819
	// Get PETSc's mat row index.
Thomas Witkowski's avatar
Thomas Witkowski committed
820
	int rowIndex = interiorMap->getMatIndex(rowComp, globalRowDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
821 822 823 824 825 826 827 828

	cols.clear();
	values.clear();

	for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); 
	     icursor != icend; ++icursor) {

	  // Global index of the current column index.
829
	  MultiIndex colMultiIndex;
Thomas Witkowski's avatar
Thomas Witkowski committed
830
	  if ((*interiorMap)[colComp].find(col(*icursor), colMultiIndex) == false)
831 832 833
	    continue;

	  int globalColDof = colMultiIndex.global;
Thomas Witkowski's avatar
Thomas Witkowski committed
834
	  // Test if the current col dof is a periodic dof.
835
	  bool periodicCol = perMap.isPeriodic(colFe, globalColDof);
836
	  // Get PETSc's mat col index.
Thomas Witkowski's avatar
Thomas Witkowski committed
837
	  int colIndex = interiorMap->getMatIndex(colComp, globalColDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
838 839 840 841 842 843 844

	  // Ignore all zero entries, expect it is a diagonal entry.
 	  if (value(*icursor) == 0.0 && rowIndex != colIndex)
 	    continue;

	  if (!periodicCol) {
	    // Calculate the exact position of the column index in the PETSc matrix.
Thomas Witkowski's avatar
Thomas Witkowski committed
845 846
 	    cols.push_back(colIndex);
 	    values.push_back(value(*icursor));
Thomas Witkowski's avatar
Thomas Witkowski committed
847 848 849 850 851
	  } else {
	    // === Row index is not periodic, but column index is. ===

	    // Create set of all periodic associations of the column index.
	    std::set<int> perAsc;
852 853 854
	    perMap.fillAssociations(colFe, globalColDof, 
				    meshDistributor->getElementObjectDb(), perAsc);

Thomas Witkowski's avatar
Thomas Witkowski committed
855 856 857 858 859 860 861 862 863
	    // Scale value to the number of periodic associations of the column index.
	    double scaledValue = 
	      value(*icursor) * pow(0.5, static_cast<double>(perAsc.size()));

	    
	    // === Create set of all matrix column indices due to the periodic ===
	    // === associations of the column DOF index.                       ===

	    vector<int> newCols;
864
	    perMap.mapDof(colFe, globalColDof, perAsc, newCols);
Thomas Witkowski's avatar
Thomas Witkowski committed
865
	    for (unsigned int i = 0; i < newCols.size(); i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
866
	      cols.push_back(interiorMap->getMatIndex(colComp, newCols[i]));
Thomas Witkowski's avatar
Thomas Witkowski committed
867 868 869 870 871
	      values.push_back(scaledValue);	      
	    }
	  }
	}

872
  	MatSetValues(getMatInterior(), 1, &rowIndex, cols.size(), 
Thomas Witkowski's avatar
Thomas Witkowski committed
873
  		     &(cols[0]), &(values[0]), ADD_VALUES);	
Thomas Witkowski's avatar
Thomas Witkowski committed
874 875 876 877 878 879 880 881 882 883 884 885 886 887
      } else {
	// === Row DOF index is periodic. ===

	// Because this row is periodic, we will have to add the entries of this 
	// matrix row to multiple rows. The following maps store to each row an
	// array of column indices and values of the entries that must be added to
	// the PETSc matrix.
	map<int, vector<int> > colsMap;
	map<int, vector<double> > valsMap;

	// Traverse all column entries.
	for (icursor_type icursor = begin<nz>(cursor), icend = end<nz>(cursor); 
	     icursor != icend; ++icursor) {
	  // Global index of the current column index.
Thomas Witkowski's avatar
Thomas Witkowski committed
888
	  int globalColDof = (*interiorMap)[colComp][col(*icursor)].global;
Thomas Witkowski's avatar
Thomas Witkowski committed
889 890 891 892 893 894 895 896 897

	  // Ignore all zero entries, expect it is a diagonal entry.
 	  if (value(*icursor) == 0.0 && globalRowDof != globalColDof)
 	    continue;

	  // === Add all periodic associations of both, the row and the column ===
	  // === indices to the set perAsc.                                    ===

	  std::set<int> perAsc;
898 899 900 901
	  perMap.fillAssociations(colFe, globalColDof, 
				  meshDistributor->getElementObjectDb(), perAsc);
	  perMap.fillAssociations(rowFe, globalRowDof, 
				  meshDistributor->getElementObjectDb(), perAsc);
Thomas Witkowski's avatar
Thomas Witkowski committed
902 903 904 905 906 907 908 909 910 911

	  // Scale the value with respect to the number of periodic associations.
	  double scaledValue = 
	    value(*icursor) * pow(0.5, static_cast<double>(perAsc.size()));


	  // === Create all matrix entries with respect to the periodic  ===
	  // === associations of the row and column indices.             ===

	  vector<pair<int, int> > entry;
Thomas Witkowski's avatar
Blbu  
Thomas Witkowski committed
912 913
	  perMap.mapDof(rowFe, colFe, make_pair(globalRowDof, globalColDof),
			perAsc, entry);
Thomas Witkowski's avatar
Thomas Witkowski committed
914 915 916

	  // === Translate the matrix entries to PETSc's matrix.

917
	  for (unsigned int i = 0; i < entry.size(); i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
918 919
	    int rowIdx = interiorMap->getMatIndex(rowComp, entry[i].first);
	    int colIdx = interiorMap->getMatIndex(colComp, entry[i].second);
Thomas Witkowski's avatar
Thomas Witkowski committed
920

921 922
	    colsMap[rowIdx].push_back(colIdx);
	    valsMap[rowIdx].push_back(scaledValue);
Thomas Witkowski's avatar
Thomas Witkowski committed
923 924 925 926 927 928 929 930 931 932 933 934
	  }
	}


	// === Finally, add all periodic rows to the PETSc matrix. ===

	for (map<int, vector<int> >::iterator rowIt = colsMap.begin();
	     rowIt != colsMap.end(); ++rowIt) {
	  TEST_EXIT_DBG(rowIt->second.size() == valsMap[rowIt->first].size())
	    ("Should not happen!\n");

	  int rowIndex = rowIt->first;
935

936
	  MatSetValues(getMatInterior(), 1, &rowIndex, rowIt->second.size(),
Thomas Witkowski's avatar
Thomas Witkowski committed
937 938 939 940 941 942 943
		       &(rowIt->second[0]), &(valsMap[rowIt->first][0]), ADD_VALUES);
	}
      }
    }
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
944 945
  void PetscSolverGlobalMatrix::setDofVector(Vec vecInterior, 
					     Vec vecCoarse,
946
					     DOFVector<double>* vec, 
Thomas Witkowski's avatar
Thomas Witkowski committed
947
					     int rowComp, 
948
					     bool rankOnly)
Thomas Witkowski's avatar
Thomas Witkowski committed
949 950 951
  {
    FUNCNAME("PetscSolverGlobalMatrix::setDofVector()");

952
    const FiniteElemSpace *feSpace = vec->getFeSpace();
953
    PeriodicMap &perMap = meshDistributor->getPeriodicMap();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
954

955
    ParallelDofMapping *rowCoarseSpace = 
Thomas Witkowski's avatar
Thomas Witkowski committed
956 957 958
      (coarseSpaceMap.size() ? coarseSpaceMap[rowComp] : NULL);

    map<DegreeOfFreedom, double> &dirichletValues = vec->getDirichletValues();
959

Thomas Witkowski's avatar
Thomas Witkowski committed
960 961 962
    // Traverse all used DOFs in the dof vector.
    DOFVector<double>::Iterator dofIt(vec, USED_DOFS);
    for (dofIt.reset(); !dofIt.end(); ++dofIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
963 964 965 966 967
      
      DegreeOfFreedom dof = dofIt.getDOFIndex();
      
      if (rankOnly && !(*interiorMap)[rowComp].isRankDof(dof))
	continue;
968

Thomas Witkowski's avatar
Thomas Witkowski committed
969 970
      // Dirichlet rows can be set only be the owner ranks.
      if (dirichletValues.count(dof) && !((*interiorMap)[rowComp].isRankDof(dof)))
Thomas Witkowski's avatar
Thomas Witkowski committed
971 972
	continue;

Thomas Witkowski's avatar
Thomas Witkowski committed
973
      if (isCoarseSpace(rowComp, dof)) {
974 975
	TEST_EXIT_DBG(vecCoarse != PETSC_NULL)("Should not happen!\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
976
	int index = rowCoarseSpace->getMatIndex(rowComp, dof);