PetscSolverGlobalMatrix.cc 22.7 KB
Newer Older
Thomas Witkowski's avatar
Thomas Witkowski committed
1 2 3 4 5 6 7 8 9 10 11 12
//
// 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.


13
#include "AMDiS.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 {

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
20
  void PetscSolverGlobalMatrix::fillPetscMatrix(Matrix<DOFMatrix*> *seqMat)
Thomas Witkowski's avatar
Thomas Witkowski committed
21 22 23
  {
    FUNCNAME("PetscSolverGlobalMatrix::fillPetscMatrix()");

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
24 25 26 27 28 29
    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();

30
    createMatVec(*seqMat);
31

Thomas Witkowski's avatar
Thomas Witkowski committed
32
    if (coarseSpaceMap.size()) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
33
      fillPetscMatrixWithCoarseSpace(seqMat);
34 35 36
      return;
    }

37 38
    // === Create PETSc vector (solution and a temporary vector). ===

Thomas Witkowski's avatar
Thomas Witkowski committed
39 40 41 42 43 44
#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
45
    int nComponents = seqMat->getNumRows();
Thomas Witkowski's avatar
Thomas Witkowski committed
46 47
    for (int i = 0; i < nComponents; i++)
      for (int j = 0; j < nComponents; j++)
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
48 49
	if ((*seqMat)[i][j])
	  setDofMatrix((*seqMat)[i][j], i, j);
Thomas Witkowski's avatar
Thomas Witkowski committed
50 51 52 53 54

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

55
    matAssembly();
56 57 58
    
    if (printMatInfo) {
      MatInfo matInfo;
59
      MatGetInfo(getMatInterior(), MAT_GLOBAL_SUM, &matInfo);
60 61 62 63 64 65 66
      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));
    }
67

68
    // === Init PETSc solver. ===
69

70 71
    KSPCreate(mpiCommGlobal, &kspInterior);
    KSPGetPC(kspInterior, &pcInterior);
72
    KSPSetOperators(kspInterior, getMatInterior(), getMatInterior(), 
73
		    SAME_NONZERO_PATTERN); 
74 75 76 77
    KSPSetTolerances(kspInterior, 0.0, 1e-8, PETSC_DEFAULT, PETSC_DEFAULT);
    KSPSetType(kspInterior, KSPBCGS);
    KSPSetOptionsPrefix(kspInterior, kspPrefix.c_str());
    KSPSetFromOptions(kspInterior);
78

79
    initPreconditioner(pcInterior);
80

81 82
    // Do not delete the solution vector, use it for the initial guess.
    if (!zeroStartVector)
83
      KSPSetInitialGuessNonzero(kspInterior, PETSC_TRUE);
84

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


91 92 93 94 95 96 97 98
  void PetscSolverGlobalMatrix::fillPetscMatrix(DOFMatrix *mat)
  {
    Matrix<DOFMatrix*> m(1, 1);
    m[0][0] = mat;
    fillPetscMatrix(&m);
  }


Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
99
  void PetscSolverGlobalMatrix::fillPetscMatrixWithCoarseSpace(Matrix<DOFMatrix*> *seqMat)
100 101
  {
    FUNCNAME("PetscSolverGlobalMatrix::fillPetscMatrixWithCoarseSpace()");
102 103

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

107
    vector<const FiniteElemSpace*> feSpaces = AMDiS::getComponentFeSpaces(*seqMat);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
108

109 110
    int nRowsRankInterior = interiorMap->getRankDofs();
    int nRowsOverallInterior = interiorMap->getOverallDofs();
111

112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130
    // === 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
131
    int nComponents = seqMat->getSize();
Thomas Witkowski's avatar
Thomas Witkowski committed
132 133 134 135 136
    for (int rowComponent = 0; rowComponent < nComponents; rowComponent++) {
      for (int colComponent = 0; colComponent < nComponents; colComponent++) {
	DOFMatrix* dofMat = (*seqMat)[rowComponent][colComponent];

	if (!dofMat)
137 138
	  continue;

Thomas Witkowski's avatar
Thomas Witkowski committed
139 140
	ParallelDofMapping *rowCoarseSpace = coarseSpaceMap[rowComponent];
	ParallelDofMapping *colCoarseSpace = coarseSpaceMap[colComponent];
Thomas Witkowski's avatar
Thomas Witkowski committed
141

Thomas Witkowski's avatar
Thomas Witkowski committed
142 143
	traits::col<Matrix>::type col(dofMat->getBaseMatrix());
	traits::const_value<Matrix>::type value(dofMat->getBaseMatrix());
144 145
	
	// Traverse all rows.
Thomas Witkowski's avatar
Thomas Witkowski committed
146 147
	for (cursor_type cursor = begin<row>(dofMat->getBaseMatrix()), 
	       cend = end<row>(dofMat->getBaseMatrix()); cursor != cend; ++cursor) {
148

Thomas Witkowski's avatar
Thomas Witkowski committed
149
	  bool isRowCoarse = 
Thomas Witkowski's avatar
Thomas Witkowski committed
150
	    isCoarseSpace(rowComponent, feSpaces[rowComponent], *cursor);
151 152 153 154 155 156 157 158 159 160
  
	  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
161 162
	    bool isColCoarse = 
	      isCoarseSpace(colComponent, feSpaces[colComponent], col(*icursor));
163

Thomas Witkowski's avatar
Thomas Witkowski committed
164 165
	    if (isColCoarse) {
	      if (isRowCoarse) {
166 167 168 169 170 171 172
		cols.push_back(col(*icursor));
		values.push_back(value(*icursor));
	      } else {
		colsOther.push_back(col(*icursor));
		valuesOther.push_back(value(*icursor));
	      }
	    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
173
	      if (isRowCoarse) {
174 175 176 177 178 179 180 181 182 183 184 185
		colsOther.push_back(col(*icursor));
		valuesOther.push_back(value(*icursor));
	      } else {
		cols.push_back(col(*icursor));
		values.push_back(value(*icursor));
	      }
	    }
	  }  // for each nnz in row


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

Thomas Witkowski's avatar
Thomas Witkowski committed
186
	  if (isRowCoarse) {
Thomas Witkowski's avatar
Thomas Witkowski committed
187 188 189
	    int rowIndex = rowCoarseSpace->getMatIndex(rowComponent, *cursor);
	    for (unsigned int i = 0; i < cols.size(); i++)
	      cols[i] = colCoarseSpace->getMatIndex(colComponent, cols[i]);
190

Thomas Witkowski's avatar
Thomas Witkowski committed
191
	    MatSetValues(getMatCoarseByComponent(rowComponent, colComponent),
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
192
			 1, &rowIndex, cols.size(),
193 194 195 196
			 &(cols[0]), &(values[0]), ADD_VALUES);

	    if (colsOther.size()) {
	      if (subdomainLevel == 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
197 198 199
		for (unsigned int i = 0; i < colsOther.size(); i++)
		  colsOther[i] = 
		    interiorMap->getMatIndex(colComponent, colsOther[i]);
200
	      } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
201 202 203 204
		for (unsigned int i = 0; i < colsOther.size(); i++)
		  colsOther[i] = 
		    interiorMap->getMatIndex(colComponent, colsOther[i]) + 
		    rStartInterior;
205
	      }
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 212
 			   &(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
	    }
	  } else {
	    int localRowIndex = 
Thomas Witkowski's avatar
Thomas Witkowski committed
213 214 215
	      (subdomainLevel == 0 ? 
	       interiorMap->getLocalMatIndex(rowComponent, *cursor) :
	       interiorMap->getMatIndex(rowComponent, *cursor));
216

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

	      if (subdomainLevel != 0)
		globalRowIndex += rStartInterior;
Thomas Witkowski's avatar
Thomas Witkowski committed
232
	      
Thomas Witkowski's avatar
Thomas Witkowski committed
233 234 235
	      for (unsigned int i = 0; i < colsOther.size(); i++)
		colsOther[i] = 
		  colCoarseSpace->getMatIndex(colComponent, colsOther[i]);
236

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

246
    matAssembly();
247

248 249 250
    // === Create solver for the non primal (thus local) variables. ===

    KSPCreate(mpiCommLocal, &kspInterior);
251
    KSPSetOperators(kspInterior, getMatInterior(), getMatInterior(),
252
		    SAME_NONZERO_PATTERN);
253 254 255 256 257 258 259 260 261 262 263 264
    KSPSetOptionsPrefix(kspInterior, "interior_");
    KSPSetType(kspInterior, KSPPREONLY);
    KSPGetPC(kspInterior, &pcInterior);
    PCSetType(pcInterior, PCLU);
    if (subdomainLevel == 0)
      PCFactorSetMatSolverPackage(pcInterior, MATSOLVERUMFPACK);
    else
      PCFactorSetMatSolverPackage(pcInterior, MATSOLVERMUMPS);
    KSPSetFromOptions(kspInterior);  
  }


265 266 267
  void PetscSolverGlobalMatrix::fillPetscRhs(SystemVector *vec)
  {
    FUNCNAME("PetscSolverGlobalMatrix::fillPetscRhs()");
Thomas Witkowski's avatar
Thomas Witkowski committed
268

269
    TEST_EXIT_DBG(vec)("No DOF vector defined!\n");
270
    TEST_EXIT_DBG(interiorMap)("No parallel DOF map defined!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
271
    
Thomas Witkowski's avatar
Thomas Witkowski committed
272
    // === Transfer values from DOF vector to the PETSc vector. === 
Thomas Witkowski's avatar
Thomas Witkowski committed
273
    if (coarseSpaceMap.size()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
274
      for (int i = 0; i < vec->getSize(); i++)
275 276
	setDofVector(getVecRhsInterior(), 
		     getVecSolCoarseByComponent(i), vec->getDOFVector(i), i);
Thomas Witkowski's avatar
Thomas Witkowski committed
277 278
    } else {
      for (int i = 0; i < vec->getSize(); i++)
279
	setDofVector(getVecRhsInterior(), vec->getDOFVector(i), i);
Thomas Witkowski's avatar
Thomas Witkowski committed
280 281
    }

282
    vecRhsAssembly();
Thomas Witkowski's avatar
Thomas Witkowski committed
283 284 285
  }


286 287
  void PetscSolverGlobalMatrix::solvePetscMatrix(SystemVector &vec, 
						 AdaptInfo *adaptInfo)
Thomas Witkowski's avatar
Thomas Witkowski committed
288 289 290 291 292 293 294
  {
    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
295 296
      TEST_EXIT(coarseSpaceMap.size() == 0)("Not yet supported!\n");

297
      VecSet(getVecSolInterior(), 0.0);
Thomas Witkowski's avatar
Thomas Witkowski committed
298 299
      
      for (int i = 0; i < nComponents; i++)
300
	setDofVector(getVecSolInterior(), vec.getDOFVector(i), i, true);
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
301

302
      vecSolAssembly();
Thomas Witkowski's avatar
Thomas Witkowski committed
303 304
    }

305

Thomas Witkowski's avatar
Thomas Witkowski committed
306 307
    MatNullSpace matNullspace;
    Vec nullspaceBasis;
308 309 310
    if (nullspace.size() > 0 || 
	hasConstantNullspace ||
	constNullspaceComponent.size() > 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
311
      TEST_EXIT_DBG(nullspace.size() <= 1)("Not yet implemented!\n");
312

313 314 315 316 317 318 319 320 321 322
      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
323
      if (nullspace.size() > 0) {
324
	VecDuplicate(getVecSolInterior(), &nullspaceBasis);
Thomas Witkowski's avatar
Thomas Witkowski committed
325 326 327 328 329
	for (int i = 0; i < nComponents; i++)
	  setDofVector(nullspaceBasis, nullspace[0]->getDOFVector(i), i, true);
	
	VecAssemblyBegin(nullspaceBasis);
	VecAssemblyEnd(nullspaceBasis);
330 331

	VecNormalize(nullspaceBasis, PETSC_NULL);
Thomas Witkowski's avatar
Thomas Witkowski committed
332 333 334 335
	
	MatNullSpaceCreate(mpiCommGlobal, (hasConstantNullspace ? PETSC_TRUE : PETSC_FALSE), 
			   1, &nullspaceBasis, &matNullspace);

336
	MatMult(getMatInterior(), nullspaceBasis, getVecSolInterior());
Thomas Witkowski's avatar
Thomas Witkowski committed
337
	PetscReal n;
338
	VecNorm(getVecSolInterior(), NORM_2, &n);
Thomas Witkowski's avatar
Thomas Witkowski committed
339 340 341 342
	MSG("NORM IS: %e\n", n);
      } else {
	MatNullSpaceCreate(mpiCommGlobal, PETSC_TRUE, 0, PETSC_NULL, &matNullspace);
      }
343

344
      MatSetNullSpace(getMatInterior(), matNullspace);
Thomas Witkowski's avatar
Thomas Witkowski committed
345
      KSPSetNullSpace(kspInterior, matNullspace);
346

Thomas Witkowski's avatar
Thomas Witkowski committed
347
      // === Remove null space, if requested. ===
348

Thomas Witkowski's avatar
Thomas Witkowski committed
349
      if (removeRhsNullspace) {
Thomas Witkowski's avatar
Thomas Witkowski committed
350
	TEST_EXIT_DBG(coarseSpaceMap.empty())("Not supported!\n");
351 352

	MSG("Remove nullspace from rhs vector.\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
353
	
354
	MatNullSpaceRemove(matNullspace, getVecRhsInterior(), PETSC_NULL);
Thomas Witkowski's avatar
Thomas Witkowski committed
355 356 357 358
      }
    } else {
      TEST_EXIT(removeRhsNullspace == false)
	("No nullspace provided that should be removed from rhs!\n");
359 360
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
361
    // PETSc.
362
    solve(getVecRhsInterior(), getVecSolInterior());
Thomas Witkowski's avatar
Thomas Witkowski committed
363

364 365

    if (nullspace.size() > 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
366
      MatNullSpaceDestroy(&matNullspace);
367 368 369 370
      VecDestroy(&nullspaceBasis);
    }


Thomas Witkowski's avatar
Thomas Witkowski committed
371 372
    // === Transfere values from PETSc's solution vectors to the DOF vectors. ===
    PetscScalar *vecPointer;
373
    VecGetArray(getVecSolInterior(), &vecPointer);    
Thomas Witkowski's avatar
Thomas Witkowski committed
374

375
    int c = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
376
    for (int i = 0; i < nComponents; i++) {
377
      DOFVector<double> &dv = *(vec.getDOFVector(i));
378

379
      DofMap& d = (*interiorMap)[dv.getFeSpace()].getMap();
380 381 382
      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
383 384
    }

385
    VecRestoreArray(getVecSolInterior(), &vecPointer);
Thomas Witkowski's avatar
Thomas Witkowski committed
386 387 388 389

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


393 394 395 396
  void PetscSolverGlobalMatrix::solveGlobal(Vec &rhs, Vec &sol)
  {
    FUNCNAME("PetscSolverGlobalMatrix::solveGlobal()");

397 398 399
    double wtime = MPI::Wtime();
    double t0 = 0.0, t1 = 0.0;

400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418
    Vec tmp;
    if (mpiCommLocal.Get_size() == 1)
      VecCreateSeq(mpiCommLocal, interiorMap->getRankDofs(), &tmp);
    else
      VecCreateMPI(mpiCommLocal,
		   interiorMap->getRankDofs(),
		   interiorMap->getOverallDofs(),
		   &tmp);

    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);

419 420 421
    t0 = MPI::Wtime() - wtime;

    wtime = MPI::Wtime();
422
    KSPSolve(kspInterior, tmp, tmp);
423 424 425
    t1 = MPI::Wtime() - wtime;

    wtime = MPI::Wtime();
426 427 428 429 430 431 432 433 434 435 436

    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);
437 438
    t0 += MPI::Wtime() - wtime;

439
    //    MSG("TIMEING: %.5f %.5f\n", t0, t1);
440 441 442
  }


443 444 445 446
  void PetscSolverGlobalMatrix::destroyMatrixData()
  {
    FUNCNAME("PetscSolverGlobalMatrix::destroyMatrixData()");

447 448
    exitPreconditioner(pcInterior);

449
    matDestroy();
Thomas Witkowski's avatar
Thomas Witkowski committed
450

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
451
    KSPDestroy(&kspInterior);
Thomas Witkowski's avatar
Thomas Witkowski committed
452 453 454
  }


455 456 457 458
  void PetscSolverGlobalMatrix::destroyVectorData()
  {
    FUNCNAME("PetscSolverGlobalMatrix::destroyVectorData()");

459
    vecDestroy();
460 461 462
  }


463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481
  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());

482
      TEST_EXIT(nComponents > 0)("No IS block for block %d defined!\n", i);
483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498

      // 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());
      }

      IS is;
      interiorMap->createIndexSet(is, blockComponents[0], nComponents);
      PCFieldSplitSetIS(pc, isNames[i].c_str(), is);
      ISDestroy(&is);
    }
  }


499 500 501 502 503 504 505 506 507 508 509 510 511 512 513
  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
514
  void PetscSolverGlobalMatrix::setDofMatrix(DOFMatrix* seqMat,
515
					     int nRowMat, int nColMat)
Thomas Witkowski's avatar
Thomas Witkowski committed
516 517 518
  {
    FUNCNAME("PetscSolverGlobalMatrix::setDofMatrix()");

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
519
    TEST_EXIT(seqMat)("No DOFMatrix!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
520 521

    using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
522
    namespace traits = mtl::traits;
Thomas Witkowski's avatar
Thomas Witkowski committed
523 524
    typedef DOFMatrix::base_matrix_type Matrix;

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
525 526
    traits::col<Matrix>::type col(seqMat->getBaseMatrix());
    traits::const_value<Matrix>::type value(seqMat->getBaseMatrix());
Thomas Witkowski's avatar
Thomas Witkowski committed
527 528 529 530 531 532 533 534 535 536 537

    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;

538 539 540
    // Get periodic mapping object
    PeriodicMap &perMap = meshDistributor->getPeriodicMap();

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
541 542
    const FiniteElemSpace *rowFe = seqMat->getRowFeSpace();
    const FiniteElemSpace *colFe = seqMat->getColFeSpace();
Thomas Witkowski's avatar
Thomas Witkowski committed
543 544 545 546
    DofMap& rowMap = (*interiorMap)[rowFe].getMap();
    DofMap& colMap = (*interiorMap)[colFe].getMap();

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

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
549 550
    for (cursor_type cursor = begin<row>(seqMat->getBaseMatrix()), 
	   cend = end<row>(seqMat->getBaseMatrix()); cursor != cend; ++cursor) {
Thomas Witkowski's avatar
Thomas Witkowski committed
551
      // Global index of the current row DOF.
Thomas Witkowski's avatar
Thomas Witkowski committed
552
      int globalRowDof = rowMap[*cursor].global;
553

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

Thomas Witkowski's avatar
Thomas Witkowski committed
557 558 559
      if (!periodicRow) {
	// === Row DOF index is not periodic. ===

560
	// Get PETSc's mat row index.
561
	int rowIndex = interiorMap->getMatIndex(nRowMat, globalRowDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
562 563 564 565 566 567 568 569

	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.
Thomas Witkowski's avatar
Thomas Witkowski committed
570
	  int globalColDof = colMap[col(*icursor)].global;
Thomas Witkowski's avatar
Thomas Witkowski committed
571
	  // Test if the current col dof is a periodic dof.
572
	  bool periodicCol = perMap.isPeriodic(colFe, globalColDof);
573
	  // Get PETSc's mat col index.
574
	  int colIndex = interiorMap->getMatIndex(nColMat, globalColDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
575 576 577 578 579 580 581

	  // 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
582 583
 	    cols.push_back(colIndex);
 	    values.push_back(value(*icursor));
Thomas Witkowski's avatar
Thomas Witkowski committed
584 585 586 587 588
	  } else {
	    // === Row index is not periodic, but column index is. ===

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

Thomas Witkowski's avatar
Thomas Witkowski committed
592 593 594 595 596 597 598 599 600
	    // 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;
601
	    perMap.mapDof(colFe, globalColDof, perAsc, newCols);
Thomas Witkowski's avatar
Thomas Witkowski committed
602
	    for (unsigned int i = 0; i < newCols.size(); i++) {
603
	      cols.push_back(interiorMap->getMatIndex(nColMat, newCols[i]));
Thomas Witkowski's avatar
Thomas Witkowski committed
604 605 606 607 608
	      values.push_back(scaledValue);	      
	    }
	  }
	}

609
  	MatSetValues(getMatInterior(), 1, &rowIndex, cols.size(), 
Thomas Witkowski's avatar
Thomas Witkowski committed
610
  		     &(cols[0]), &(values[0]), ADD_VALUES);	
Thomas Witkowski's avatar
Thomas Witkowski committed
611 612 613 614 615 616 617 618 619 620 621 622 623 624
      } 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.
625
	  int globalColDof = (*interiorMap)[colFe][col(*icursor)].global;
Thomas Witkowski's avatar
Thomas Witkowski committed
626 627 628 629 630 631 632 633 634

	  // 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;
635 636 637 638
	  perMap.fillAssociations(colFe, globalColDof, 
				  meshDistributor->getElementObjectDb(), perAsc);
	  perMap.fillAssociations(rowFe, globalRowDof, 
				  meshDistributor->getElementObjectDb(), perAsc);
Thomas Witkowski's avatar
Thomas Witkowski committed
639 640 641 642 643 644 645 646 647 648

	  // 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
649 650
	  perMap.mapDof(rowFe, colFe, make_pair(globalRowDof, globalColDof),
			perAsc, entry);
Thomas Witkowski's avatar
Thomas Witkowski committed
651 652 653

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

654
	  for (unsigned int i = 0; i < entry.size(); i++) {
655 656
	    int rowIdx = interiorMap->getMatIndex(nRowMat, entry[i].first);
	    int colIdx = interiorMap->getMatIndex(nColMat, entry[i].second);
Thomas Witkowski's avatar
Thomas Witkowski committed
657

658 659
	    colsMap[rowIdx].push_back(colIdx);
	    valsMap[rowIdx].push_back(scaledValue);
Thomas Witkowski's avatar
Thomas Witkowski committed
660 661 662 663 664 665 666 667 668 669 670 671
	  }
	}


	// === 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;
672

673
	  MatSetValues(getMatInterior(), 1, &rowIndex, rowIt->second.size(),
Thomas Witkowski's avatar
Thomas Witkowski committed
674 675 676 677 678 679 680
		       &(rowIt->second[0]), &(valsMap[rowIt->first][0]), ADD_VALUES);
	}
      }
    }
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
681 682
  void PetscSolverGlobalMatrix::setDofVector(Vec vecInterior, 
					     Vec vecCoarse,
683
					     DOFVector<double>* vec, 
684
					     int nRowVec, 
685
					     bool rankOnly)
Thomas Witkowski's avatar
Thomas Witkowski committed
686 687 688
  {
    FUNCNAME("PetscSolverGlobalMatrix::setDofVector()");

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

692 693
    ParallelDofMapping *rowCoarseSpace = 
      (coarseSpaceMap.size() ? coarseSpaceMap[nRowVec] : NULL);
694

Thomas Witkowski's avatar
Thomas Witkowski committed
695 696 697
    // Traverse all used DOFs in the dof vector.
    DOFVector<double>::Iterator dofIt(vec, USED_DOFS);
    for (dofIt.reset(); !dofIt.end(); ++dofIt) {
698

699
      if (rankOnly && !(*interiorMap)[feSpace].isRankDof(dofIt.getDOFIndex()))
Thomas Witkowski's avatar
Thomas Witkowski committed
700 701
	continue;

Thomas Witkowski's avatar
Thomas Witkowski committed
702
      if (isCoarseSpace(nRowVec, feSpace, dofIt.getDOFIndex())) {
703 704
	TEST_EXIT_DBG(vecCoarse != PETSC_NULL)("Should not happen!\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
705
	int index = rowCoarseSpace->getMatIndex(nRowVec, dofIt.getDOFIndex());
706 707 708 709 710 711 712 713 714 715 716 717 718 719
	VecSetValue(vecCoarse, index, *dofIt, ADD_VALUES);
      } else {
	// Calculate global row index of the DOF.
	DegreeOfFreedom globalRowDof = 
	  (*interiorMap)[feSpace][dofIt.getDOFIndex()].global;
	
	// Get PETSc's mat index of the row DOF.
	int index = 0;
	if (interiorMap->isMatIndexFromGlobal())
	  index = 
	    interiorMap->getMatIndex(nRowVec, globalRowDof) + rStartInterior;
	else
	  index =
	    interiorMap->getMatIndex(nRowVec, dofIt.getDOFIndex()) + rStartInterior;
720

721
	if (perMap.isPeriodic(feSpace, globalRowDof)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
722 723 724
	  std::set<int>& perAsc = perMap.getAssociations(feSpace, globalRowDof);
	  double value = *dofIt / (perAsc.size() + 1.0);
	  VecSetValue(vecInterior, index, value, ADD_VALUES);
725

Thomas Witkowski's avatar
Thomas Witkowski committed
726 727 728 729
	  for (std::set<int>::iterator perIt = perAsc.begin(); 
	       perIt != perAsc.end(); ++perIt) {
	    int mappedDof = perMap.map(feSpace, *perIt, globalRowDof);
	    int mappedIndex = interiorMap->getMatIndex(nRowVec, mappedDof);
730

Thomas Witkowski's avatar
Thomas Witkowski committed
731
	    VecSetValue(vecInterior, mappedIndex, value, ADD_VALUES);
732 733 734 735
	  }	  
	} else {	  
	  // The DOF index is not periodic.
	  
Thomas Witkowski's avatar
Thomas Witkowski committed
736 737
	  VecSetValue(vecInterior, index, *dofIt, ADD_VALUES);
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
738 739 740 741 742
      }
    }
  }

}