PetscSolverGlobalMatrix.cc 30.1 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 {

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

24
    if (coarseSpaceMap != NULL) {
Thomas Witkowski's avatar
Thomas Witkowski committed
25
      updateSubdomainData();
26
27
28
29
      fillPetscMatrixWithCoarseSpace(mat);
      return;
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
30
    TEST_EXIT_DBG(meshDistributor)("No mesh distributor object defined!\n");
31
    TEST_EXIT_DBG(interiorMap)("No parallel mapping object defined!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
32
    TEST_EXIT_DBG(mat)("No DOF matrix defined!\n");
33
    
Thomas Witkowski's avatar
Thomas Witkowski committed
34
    double wtime = MPI::Wtime();
35

36
    // === If required, recompute non zero structure of the matrix. ===
Thomas Witkowski's avatar
Thomas Witkowski committed
37

38
    if (checkMeshChange(mat))
Thomas Witkowski's avatar
Thomas Witkowski committed
39
      nnzInterior.create(mat, mpiCommGlobal, *interiorMap, 
40
41
			 &(meshDistributor->getPeriodicMap()),
			 meshDistributor->getElementObjectDb());
Thomas Witkowski's avatar
Thomas Witkowski committed
42

43
44
    // === Create PETSc vector (solution and a temporary vector). ===

45
46
    int nRankRows = interiorMap->getRankDofs();
    int nOverallRows = interiorMap->getOverallDofs();
47

48
    VecCreateMPI(mpiCommGlobal, nRankRows, nOverallRows, &petscSolVec);
49
50


Thomas Witkowski's avatar
Thomas Witkowski committed
51
52
    // === Create PETSc matrix with the computed nnz data structure. ===

53
    MatCreateAIJ(mpiCommGlobal, nRankRows, nRankRows, 
54
		 nOverallRows, nOverallRows,
55
56
57
		 0, nnzInterior.dnnz, 
		 0, nnzInterior.onnz, 
		 &matIntInt);
58
59

    MatSetOption(matIntInt, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
60

Thomas Witkowski's avatar
Thomas Witkowski committed
61
62
63
64
65
66
#if (DEBUG != 0)
    MSG("Fill petsc matrix 1 needed %.5f seconds\n", MPI::Wtime() - wtime);
#endif

#if (DEBUG != 0)
    int a, b;
67
68
69
    MatGetOwnershipRange(matIntInt, &a, &b);
    TEST_EXIT(a == interiorMap->getStartDofs())("Wrong matrix ownership range!\n");
    TEST_EXIT(b == interiorMap->getStartDofs() + nRankRows)
Thomas Witkowski's avatar
Thomas Witkowski committed
70
71
72
73
74
75
      ("Wrong matrix ownership range!\n");
#endif


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

76
    int nComponents = mat->getNumRows();
Thomas Witkowski's avatar
Thomas Witkowski committed
77
78
79
    for (int i = 0; i < nComponents; i++)
      for (int j = 0; j < nComponents; j++)
	if ((*mat)[i][j])
80
	  setDofMatrix((*mat)[i][j], i, j);
Thomas Witkowski's avatar
Thomas Witkowski committed
81
82
83
84
85

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

86
87
    MatAssemblyBegin(matIntInt, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(matIntInt, MAT_FINAL_ASSEMBLY);
Thomas Witkowski's avatar
Thomas Witkowski committed
88

89
90
91
92
93
94
95
96
97
98
99
    
    if (printMatInfo) {
      MatInfo matInfo;
      MatGetInfo(matIntInt, MAT_GLOBAL_SUM, &matInfo);
      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));
    }
100
101
102

    // === Remove Dirichlet BC DOFs. ===

103
    //    removeDirichletBcDofs(mat);
104
105
    

106
    // === Init PETSc solver. ===
107

108
109
110
111
112
113
114
    KSPCreate(mpiCommGlobal, &kspInterior);
    KSPGetPC(kspInterior, &pcInterior);
    KSPSetOperators(kspInterior, matIntInt, matIntInt, SAME_NONZERO_PATTERN); 
    KSPSetTolerances(kspInterior, 0.0, 1e-8, PETSC_DEFAULT, PETSC_DEFAULT);
    KSPSetType(kspInterior, KSPBCGS);
    KSPSetOptionsPrefix(kspInterior, kspPrefix.c_str());
    KSPSetFromOptions(kspInterior);
115

116
    initPreconditioner(pcInterior);
117

118
119
    // Do not delete the solution vector, use it for the initial guess.
    if (!zeroStartVector)
120
      KSPSetInitialGuessNonzero(kspInterior, PETSC_TRUE);
121

Thomas Witkowski's avatar
Thomas Witkowski committed
122
123
124
#if (DEBUG != 0)
    MSG("Fill petsc matrix 3 needed %.5f seconds\n", MPI::Wtime() - wtime);
#endif
125
126
127
  }


128
129
130
131
132
133
134
135
  void PetscSolverGlobalMatrix::fillPetscMatrix(DOFMatrix *mat)
  {
    Matrix<DOFMatrix*> m(1, 1);
    m[0][0] = mat;
    fillPetscMatrix(&m);
  }


136
137
138
  void PetscSolverGlobalMatrix::fillPetscMatrixWithCoarseSpace(Matrix<DOFMatrix*> *mat)
  {
    FUNCNAME("PetscSolverGlobalMatrix::fillPetscMatrixWithCoarseSpace()");
139
140

    TEST_EXIT_DBG(interiorMap)("Should not happen!\n");
141

142
143
144
145
146
    vector<const FiniteElemSpace*> feSpaces = getFeSpaces(mat);

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

147
148
149
150
    // === If required, recompute non zero structure of the matrix. ===

    bool localMatrix = (subdomainLevel == 0);
    if (checkMeshChange(mat, localMatrix)) {
151
152
153
      nnzInterior.create(mat, mpiCommGlobal, *interiorMap, NULL, 
			 meshDistributor->getElementObjectDb(),
			 localMatrix);
154
155

      if (coarseSpaceMap) {
156
157
158
	nnzCoarse.create(mat, mpiCommGlobal, *coarseSpaceMap, NULL, meshDistributor->getElementObjectDb());
	nnzCoarseInt.create(mat, mpiCommGlobal, *coarseSpaceMap, *interiorMap, NULL, meshDistributor->getElementObjectDb());
	nnzIntCoarse.create(mat, mpiCommGlobal, *interiorMap, *coarseSpaceMap, NULL, meshDistributor->getElementObjectDb());
159
      }
160
161
162
    }

    if (localMatrix) {
163
      MatCreateSeqAIJ(mpiCommLocal, nRowsRankInterior, nRowsRankInterior,
164
165
		      0, nnzInterior.dnnz, 
		      &matIntInt);
166
    } else {
167
168
169
      MatCreateAIJ(mpiCommLocal, 
		   nRowsRankInterior, nRowsRankInterior,
		   nRowsOverallInterior, nRowsOverallInterior,
170
		   0, nnzInterior.dnnz, 0, nnzInterior.onnz, 
171
		   &matIntInt);
172
    }
173
174

   
175
176
177
178
    if (coarseSpaceMap) {
      int nRowsRankCoarse = coarseSpaceMap->getRankDofs();
      int nRowsOverallCoarse = coarseSpaceMap->getOverallDofs();

179
180
181
      MatCreateAIJ(mpiCommGlobal,
		   nRowsRankCoarse, nRowsRankCoarse,
		   nRowsOverallCoarse, nRowsOverallCoarse,
182
		   0, nnzCoarse.dnnz, 0, nnzCoarse.onnz, 
183
184
		   &matCoarseCoarse);

185
186
187
      MatCreateAIJ(mpiCommGlobal,
		   nRowsRankCoarse, nRowsRankInterior,
		   nRowsOverallCoarse, nGlobalOverallInterior,
188
		   0, nnzCoarseInt.dnnz, 0, nnzCoarseInt.onnz,
189
190
		   &matCoarseInt);

191
192
193
      MatCreateAIJ(mpiCommGlobal,
		   nRowsRankInterior, nRowsRankCoarse,
		   nGlobalOverallInterior, nRowsOverallCoarse,
194
		   0, nnzIntCoarse.dnnz, 0, nnzIntCoarse.onnz,
195
		   &matIntCoarse);
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
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
    // === 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.                                       ===

    int nComponents = mat->getSize();
    for (int i = 0; i < nComponents; i++) {
      for (int j = 0; j < nComponents; j++) {
	if (!(*mat)[i][j])
	  continue;

	traits::col<Matrix>::type col((*mat)[i][j]->getBaseMatrix());
	traits::const_value<Matrix>::type value((*mat)[i][j]->getBaseMatrix());
	
	// Traverse all rows.
	for (cursor_type cursor = begin<row>((*mat)[i][j]->getBaseMatrix()), 
	       cend = end<row>((*mat)[i][j]->getBaseMatrix()); cursor != cend; ++cursor) {

	  bool rowPrimal = isCoarseSpace(feSpaces[i], *cursor);
  
	  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) {

	    bool colPrimal = isCoarseSpace(feSpaces[j], col(*icursor));

	    if (colPrimal) {
	      if (rowPrimal) {
		cols.push_back(col(*icursor));
		values.push_back(value(*icursor));
	      } else {
		colsOther.push_back(col(*icursor));
		valuesOther.push_back(value(*icursor));
	      }
	    } else {
	      if (rowPrimal) {
		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. ===

	  if (rowPrimal) {
	    int rowIndex = coarseSpaceMap->getMatIndex(i, *cursor);
	    for (unsigned int k = 0; k < cols.size(); k++)
	      cols[k] = coarseSpaceMap->getMatIndex(j, cols[k]);

	    MatSetValues(matCoarseCoarse, 1, &rowIndex, cols.size(),
			 &(cols[0]), &(values[0]), ADD_VALUES);

	    if (colsOther.size()) {
	      if (subdomainLevel == 0) {
276
		for (unsigned int k = 0; k < colsOther.size(); k++)
277
278
279
280
281
282
		  colsOther[k] = interiorMap->getMatIndex(j, colsOther[k]);
	      } else {
		for (unsigned int k = 0; k < colsOther.size(); k++)
		  colsOther[k] = 
		    interiorMap->getMatIndex(j, colsOther[k]) + rStartInterior;
	      }
283

284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
	      MatSetValues(matCoarseInt, 1, &rowIndex, colsOther.size(),
 			   &(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
	    }
	  } else {
	    int localRowIndex = 
	      (subdomainLevel == 0 ? interiorMap->getLocalMatIndex(i, *cursor) :
	       interiorMap->getMatIndex(i, *cursor));

	    for (unsigned int k = 0; k < cols.size(); k++) {
	      if (subdomainLevel == 0)
		cols[k] = interiorMap->getLocalMatIndex(j, cols[k]);
	      else
		cols[k] = interiorMap->getMatIndex(j, cols[k]);
	    }
	    
  	    MatSetValues(matIntInt, 1, &localRowIndex, cols.size(),
  			 &(cols[0]), &(values[0]), ADD_VALUES);

	    if (colsOther.size()) {
	      int globalRowIndex = interiorMap->getMatIndex(i, *cursor);

	      if (subdomainLevel != 0)
		globalRowIndex += rStartInterior;

	      for (unsigned int k = 0; k < colsOther.size(); k++)
		colsOther[k] = coarseSpaceMap->getMatIndex(j, colsOther[k]);

  	      MatSetValues(matIntCoarse, 1, &globalRowIndex, colsOther.size(),
  			   &(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
	    }
	  }
	} 
      }
    }

319

320
321
322
323
324
    // === Start global assembly procedure. ===

    MatAssemblyBegin(matIntInt, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(matIntInt, MAT_FINAL_ASSEMBLY);

325

326
327
328
329
330
331
332
333
334
335
336
    if (coarseSpaceMap) {
      MatAssemblyBegin(matCoarseCoarse, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matCoarseCoarse, MAT_FINAL_ASSEMBLY);
      
      MatAssemblyBegin(matIntCoarse, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matIntCoarse, MAT_FINAL_ASSEMBLY);
      
      MatAssemblyBegin(matCoarseInt, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(matCoarseInt, MAT_FINAL_ASSEMBLY);
    }

337

338
339
    // === Remove Dirichlet BC DOFs. ===

340
    //    removeDirichletBcDofs(mat);
341

342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
    // === Create solver for the non primal (thus local) variables. ===

    KSPCreate(mpiCommLocal, &kspInterior);
    KSPSetOperators(kspInterior, matIntInt, matIntInt, SAME_NONZERO_PATTERN);
    KSPSetOptionsPrefix(kspInterior, "interior_");
    KSPSetType(kspInterior, KSPPREONLY);
    KSPGetPC(kspInterior, &pcInterior);
    PCSetType(pcInterior, PCLU);
    if (subdomainLevel == 0)
      PCFactorSetMatSolverPackage(pcInterior, MATSOLVERUMFPACK);
    else
      PCFactorSetMatSolverPackage(pcInterior, MATSOLVERMUMPS);
    KSPSetFromOptions(kspInterior);  
  }


358
359
360
  void PetscSolverGlobalMatrix::fillPetscRhs(SystemVector *vec)
  {
    FUNCNAME("PetscSolverGlobalMatrix::fillPetscRhs()");
Thomas Witkowski's avatar
Thomas Witkowski committed
361

Thomas Witkowski's avatar
Thomas Witkowski committed
362
363
364
365
    VecCreateMPI(mpiCommGlobal, 
		 interiorMap->getRankDofs(), 
		 nGlobalOverallInterior,
		 &rhsInterior);
366

Thomas Witkowski's avatar
Thomas Witkowski committed
367
368
369
370
371
372
    if (coarseSpaceMap) 
      VecCreateMPI(mpiCommGlobal, 
		   coarseSpaceMap->getRankDofs(), 
		   coarseSpaceMap->getOverallDofs(),
		   &rhsCoarseSpace);
    
373
    TEST_EXIT_DBG(vec)("No DOF vector defined!\n");
374
    TEST_EXIT_DBG(interiorMap)("No parallel DOF map defined!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
375
    
376

Thomas Witkowski's avatar
Thomas Witkowski committed
377
    // === Transfer values from DOF vector to the PETSc vector. === 
Thomas Witkowski's avatar
Thomas Witkowski committed
378
    if (coarseSpaceMap) {
Thomas Witkowski's avatar
Thomas Witkowski committed
379
380
      for (int i = 0; i < vec->getSize(); i++)
	setDofVector(rhsInterior, rhsCoarseSpace, vec->getDOFVector(i), i);
Thomas Witkowski's avatar
Thomas Witkowski committed
381
382
383
384
    } else {
      for (int i = 0; i < vec->getSize(); i++)
	setDofVector(rhsInterior, vec->getDOFVector(i), i);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
385

386
387
    VecAssemblyBegin(rhsInterior);
    VecAssemblyEnd(rhsInterior);
388

Thomas Witkowski's avatar
Thomas Witkowski committed
389
390
391
392
393
    if (coarseSpaceMap) {
      VecAssemblyBegin(rhsCoarseSpace);
      VecAssemblyEnd(rhsCoarseSpace);
    }

394
395
396

    // === Remove Dirichlet BC DOFs. ===
    //    removeDirichletBcDofs(vec);
Thomas Witkowski's avatar
Thomas Witkowski committed
397
 
Thomas Witkowski's avatar
Thomas Witkowski committed
398
399
400
  }


401
402
  void PetscSolverGlobalMatrix::solvePetscMatrix(SystemVector &vec, 
						 AdaptInfo *adaptInfo)
Thomas Witkowski's avatar
Thomas Witkowski committed
403
404
405
406
407
408
409
410
411
412
  {
    FUNCNAME("PetscSolverGlobalMatrix::solvePetscMatrix()");

    int nComponents = vec.getSize();

    // === Set old solution to be initiual guess for PETSc solver. ===
    if (!zeroStartVector) {
      VecSet(petscSolVec, 0.0);
      
      for (int i = 0; i < nComponents; i++)
413
	setDofVector(petscSolVec, vec.getDOFVector(i), i, true);
Thomas Witkowski's avatar
Thomas Witkowski committed
414
415
416
417
418
      
      VecAssemblyBegin(petscSolVec);
      VecAssemblyEnd(petscSolVec);
    }

419

Thomas Witkowski's avatar
Thomas Witkowski committed
420
421
    MatNullSpace matNullspace;
    Vec nullspaceBasis;
422
423
424
    if (nullspace.size() > 0 || 
	hasConstantNullspace ||
	constNullspaceComponent.size() > 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
425
      TEST_EXIT_DBG(nullspace.size() <= 1)("Not yet implemented!\n");
426

427
428
429
430
431
432
433
434
435
436
      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
437
438
439
440
441
442
443
      if (nullspace.size() > 0) {
	VecDuplicate(petscSolVec, &nullspaceBasis);
	for (int i = 0; i < nComponents; i++)
	  setDofVector(nullspaceBasis, nullspace[0]->getDOFVector(i), i, true);
	
	VecAssemblyBegin(nullspaceBasis);
	VecAssemblyEnd(nullspaceBasis);
444
445

	VecNormalize(nullspaceBasis, PETSC_NULL);
Thomas Witkowski's avatar
Thomas Witkowski committed
446
447
448
449
450
451
452
453
454
455
456
	
	MatNullSpaceCreate(mpiCommGlobal, (hasConstantNullspace ? PETSC_TRUE : PETSC_FALSE), 
			   1, &nullspaceBasis, &matNullspace);

	MatMult(matIntInt, nullspaceBasis, petscSolVec);
	PetscReal n;
	VecNorm(petscSolVec, NORM_2, &n);
	MSG("NORM IS: %e\n", n);
      } else {
	MatNullSpaceCreate(mpiCommGlobal, PETSC_TRUE, 0, PETSC_NULL, &matNullspace);
      }
457

458
      MatSetNullSpace(matIntInt, matNullspace);
Thomas Witkowski's avatar
Thomas Witkowski committed
459
      KSPSetNullSpace(kspInterior, matNullspace);
460

Thomas Witkowski's avatar
Thomas Witkowski committed
461
      // === Remove null space, if requested. ===
462

Thomas Witkowski's avatar
Thomas Witkowski committed
463
464
465
466
467
468
469
470
      if (removeRhsNullspace) {
	TEST_EXIT_DBG(coarseSpaceMap == NULL)("Not supported!\n");
	
	MatNullSpaceRemove(matNullspace, rhsInterior, PETSC_NULL);
      }
    } else {
      TEST_EXIT(removeRhsNullspace == false)
	("No nullspace provided that should be removed from rhs!\n");
471
472
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
473
    // PETSc.
474
    solve(rhsInterior, petscSolVec);
Thomas Witkowski's avatar
Thomas Witkowski committed
475

476
477

    if (nullspace.size() > 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
478
      MatNullSpaceDestroy(&matNullspace);
479
480
481
482
      VecDestroy(&nullspaceBasis);
    }


Thomas Witkowski's avatar
Thomas Witkowski committed
483
484
    // === Transfere values from PETSc's solution vectors to the DOF vectors. ===
    PetscScalar *vecPointer;
485
    VecGetArray(petscSolVec, &vecPointer);    
Thomas Witkowski's avatar
Thomas Witkowski committed
486

487
    int c = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
488
    for (int i = 0; i < nComponents; i++) {
489
      DOFVector<double> &dv = *(vec.getDOFVector(i));
490

491
      DofMap& d = (*interiorMap)[dv.getFeSpace()].getMap();
492
493
494
      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
495
496
497
    }

    VecRestoreArray(petscSolVec, &vecPointer);
Thomas Witkowski's avatar
Thomas Witkowski committed
498
499
500
501

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


505
506
507
508
  void PetscSolverGlobalMatrix::solveGlobal(Vec &rhs, Vec &sol)
  {
    FUNCNAME("PetscSolverGlobalMatrix::solveGlobal()");

509
510
511
    double wtime = MPI::Wtime();
    double t0 = 0.0, t1 = 0.0;

512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
    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);

531
532
533
    t0 = MPI::Wtime() - wtime;

    wtime = MPI::Wtime();
534
    KSPSolve(kspInterior, tmp, tmp);
535
536
537
    t1 = MPI::Wtime() - wtime;

    wtime = MPI::Wtime();
538
539
540
541
542
543
544
545
546
547
548

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

551
    //    MSG("TIMEING: %.5f %.5f\n", t0, t1);
552
553
554
  }


555
556
557
558
  void PetscSolverGlobalMatrix::destroyMatrixData()
  {
    FUNCNAME("PetscSolverGlobalMatrix::destroyMatrixData()");

559
560
    exitPreconditioner(pcInterior);

561
562
    MatDestroy(&matIntInt);
    KSPDestroy(&kspInterior);
Thomas Witkowski's avatar
Thomas Witkowski committed
563

564
565
566
567
568
569
570
    if (coarseSpaceMap) {
      MatDestroy(&matCoarseCoarse);
      MatDestroy(&matCoarseInt);
      MatDestroy(&matIntCoarse);
    }

    if (petscSolVec != PETSC_NULL) {
Thomas Witkowski's avatar
Thomas Witkowski committed
571
      VecDestroy(&petscSolVec);
572
573
      petscSolVec = PETSC_NULL;
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
574
575
576
  }


577
578
579
580
  void PetscSolverGlobalMatrix::destroyVectorData()
  {
    FUNCNAME("PetscSolverGlobalMatrix::destroyVectorData()");

581
    VecDestroy(&rhsInterior);
582
583
584

    if (coarseSpaceMap)
      VecDestroy(&rhsCoarseSpace);
585
586
587
  }


588
589
590
591
592
593
594
595
596
597
  bool PetscSolverGlobalMatrix::checkMeshChange(Matrix<DOFMatrix*> *mat,
						bool localMatrix)
  {
    FUNCNAME("PetscSolverGlobalMatrix::checkMeshChange()");

    int recvAllValues = 0;
    int sendValue = 
      static_cast<int>(meshDistributor->getLastMeshChangeIndex() != lastMeshNnz);
    mpiCommGlobal.Allreduce(&sendValue, &recvAllValues, 1, MPI_INT, MPI_SUM);

598
    if (!nnzInterior.dnnz || recvAllValues != 0 || alwaysCreateNnzStructure) {
599
600
601
602
603
      vector<const FiniteElemSpace*> feSpaces = getFeSpaces(mat);
      
      interiorMap->setComputeMatIndex(true, !localMatrix);
      interiorMap->update(feSpaces);

604
      nnzInterior.clear();
605

606
607
608
609
      if (coarseSpaceMap) {
	nnzCoarse.clear();
	nnzCoarseInt.clear();
	nnzIntCoarse.clear();
610
611
612
613
614
615
616
617
618
619
620
621
      }

      updateSubdomainData();
      lastMeshNnz = meshDistributor->getLastMeshChangeIndex();

      return true;
    }

    return false;
  }


622
623
624
625
626
627
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
  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());

      TEST_EXIT(nComponents > 0)("No is block for block %d defined!\n", i);

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


658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
  void PetscSolverGlobalMatrix::initPreconditioner(PC pc)
  {
    FUNCNAME("PetscSolverGlobalMatrix::initPreconditioner()");

    PCSetFromOptions(pc);
    createFieldSplit(pc);
  }

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


673
674
  void PetscSolverGlobalMatrix::setDofMatrix(DOFMatrix* mat,
					     int nRowMat, int nColMat)
Thomas Witkowski's avatar
Thomas Witkowski committed
675
676
677
678
679
680
  {
    FUNCNAME("PetscSolverGlobalMatrix::setDofMatrix()");

    TEST_EXIT(mat)("No DOFMatrix!\n");

    using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
681
    namespace traits = mtl::traits;
Thomas Witkowski's avatar
Thomas Witkowski committed
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
    typedef DOFMatrix::base_matrix_type Matrix;

    traits::col<Matrix>::type col(mat->getBaseMatrix());
    traits::const_value<Matrix>::type value(mat->getBaseMatrix());

    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;

697
698
699
    // Get periodic mapping object
    PeriodicMap &perMap = meshDistributor->getPeriodicMap();

Thomas Witkowski's avatar
Thomas Witkowski committed
700
701
702
703
704
705
    const FiniteElemSpace *rowFe = mat->getRowFeSpace();
    const FiniteElemSpace *colFe = mat->getColFeSpace();
    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
706
707
708
709
710
    // === to the PETSc matrix.                                               ===

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
716
717
718
      if (!periodicRow) {
	// === Row DOF index is not periodic. ===

719
	// Get PETSc's mat row index.
720
	int rowIndex = interiorMap->getMatIndex(nRowMat, globalRowDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
721
722
723
724
725
726
727
728

	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
729
	  int globalColDof = colMap[col(*icursor)].global;
Thomas Witkowski's avatar
Thomas Witkowski committed
730
	  // Test if the current col dof is a periodic dof.
731
	  bool periodicCol = perMap.isPeriodic(colFe, globalColDof);
732
	  // Get PETSc's mat col index.
733
	  int colIndex = interiorMap->getMatIndex(nColMat, globalColDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
734
735
736
737
738
739
740

	  // 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
741
742
 	    cols.push_back(colIndex);
 	    values.push_back(value(*icursor));
Thomas Witkowski's avatar
Thomas Witkowski committed
743
744
745
746
747
	  } else {
	    // === Row index is not periodic, but column index is. ===

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

Thomas Witkowski's avatar
Thomas Witkowski committed
751
752
753
754
755
756
757
758
759
	    // 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;
760
	    perMap.mapDof(colFe, globalColDof, perAsc, newCols);
Thomas Witkowski's avatar
Thomas Witkowski committed
761
	    for (unsigned int i = 0; i < newCols.size(); i++) {
762
	      cols.push_back(interiorMap->getMatIndex(nColMat, newCols[i]));
Thomas Witkowski's avatar
Thomas Witkowski committed
763
764
765
766
767
	      values.push_back(scaledValue);	      
	    }
	  }
	}

Thomas Witkowski's avatar
Thomas Witkowski committed
768
769
  	MatSetValues(matIntInt, 1, &rowIndex, cols.size(), 
  		     &(cols[0]), &(values[0]), ADD_VALUES);	
Thomas Witkowski's avatar
Thomas Witkowski committed
770
771
772
773
774
775
776
777
778
779
780
781
782
783
      } 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.
784
	  int globalColDof = (*interiorMap)[colFe][col(*icursor)].global;
Thomas Witkowski's avatar
Thomas Witkowski committed
785
786
787
788
789
790
791
792
793

	  // 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;
794
795
796
797
	  perMap.fillAssociations(colFe, globalColDof, 
				  meshDistributor->getElementObjectDb(), perAsc);
	  perMap.fillAssociations(rowFe, globalRowDof, 
				  meshDistributor->getElementObjectDb(), perAsc);
Thomas Witkowski's avatar
Thomas Witkowski committed
798
799
800
801
802
803
804
805
806
807
808
809
810
811

	  // 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;
	  
	  // First, add the original entry.
	  entry.push_back(make_pair(globalRowDof, globalColDof));

812
813
	  // Then, traverse the periodic associations of the row and column
	  // indices and create the corresponding entries.
Thomas Witkowski's avatar
Thomas Witkowski committed
814
815
816
817
	  for (std::set<int>::iterator it = perAsc.begin(); it != perAsc.end(); ++it) {
	    int nEntry = static_cast<int>(entry.size());
	    for (int i = 0; i < nEntry; i++) {
	      int perRowDof = 0;
818

819
820
	      if (perMap.isPeriodic(rowFe, *it, entry[i].first))
		perRowDof = perMap.map(rowFe, *it, entry[i].first);
Thomas Witkowski's avatar
Thomas Witkowski committed
821
822
823
824
	      else
		perRowDof = entry[i].first;

	      int perColDof;
825
826
	      if (perMap.isPeriodic(colFe, *it, entry[i].second))
		perColDof = perMap.map(colFe, *it, entry[i].second);
Thomas Witkowski's avatar
Thomas Witkowski committed
827
828
829
830
831
832
833
834
835
836
	      else
		perColDof = entry[i].second;	      	      
	      
	      entry.push_back(make_pair(perRowDof, perColDof));
	    }
	  }


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

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

841
842
	    colsMap[rowIdx].push_back(colIdx);
	    valsMap[rowIdx].push_back(scaledValue);
Thomas Witkowski's avatar
Thomas Witkowski committed
843
844
845
846
847
848
849
850
851
852
853
854
	  }
	}


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

856
	  MatSetValues(matIntInt, 1, &rowIndex, rowIt->second.size(),
Thomas Witkowski's avatar
Thomas Witkowski committed
857
858
859
860
861
862
863
		       &(rowIt->second[0]), &(valsMap[rowIt->first][0]), ADD_VALUES);
	}
      }
    }
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
864
865
  void PetscSolverGlobalMatrix::setDofVector(Vec vecInterior, 
					     Vec vecCoarse,
866
					     DOFVector<double>* vec, 
867
					     int nRowVec, 
868
					     bool rankOnly)
Thomas Witkowski's avatar
Thomas Witkowski committed
869
870
871
  {
    FUNCNAME("PetscSolverGlobalMatrix::setDofVector()");

872
    const FiniteElemSpace *feSpace = vec->getFeSpace();
873
    PeriodicMap &perMap = meshDistributor->getPeriodicMap();
874

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

879
      if (rankOnly && !(*interiorMap)[feSpace].isRankDof(dofIt.getDOFIndex()))
Thomas Witkowski's avatar
Thomas Witkowski committed
880
881
	continue;

882
      if (isCoarseSpace(feSpace, dofIt.getDOFIndex())) {
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
	TEST_EXIT_DBG(vecCoarse != PETSC_NULL)("Should not happen!\n");

	int index = coarseSpaceMap->getMatIndex(nRowVec, dofIt.getDOFIndex());
	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;
900

901
	if (perMap.isPeriodic(feSpace, globalRowDof)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
902
903
904
	  std::set<int>& perAsc = perMap.getAssociations(feSpace, globalRowDof);
	  double value = *dofIt / (perAsc.size() + 1.0);
	  VecSetValue(vecInterior, index, value, ADD_VALUES);
905

Thomas Witkowski's avatar
Thomas Witkowski committed
906
907
908
909
	  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);
910

Thomas Witkowski's avatar
Thomas Witkowski committed
911
	    VecSetValue(vecInterior, mappedIndex, value, ADD_VALUES);
912
913
914
915
	  }	  
	} else {	  
	  // The DOF index is not periodic.
	  
Thomas Witkowski's avatar
Thomas Witkowski committed
916
917
	  VecSetValue(vecInterior, index, *dofIt, ADD_VALUES);
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
918
919
920
921
922
      }
    }
  }


923
924
925
926
  void PetscSolverGlobalMatrix::removeDirichletBcDofs(Matrix<DOFMatrix*> *mat)
  {
    FUNCNAME("PetscSolverGlobalMatrix::removeDirichletBcDofs()");

Thomas Witkowski's avatar
Thomas Witkowski committed
927
#if 0
928
929
930
931
    vector<int> dofsInterior, dofsCoarse;

    int nComponents = mat->getNumRows();
    for (int i = 0; i < nComponents; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
      for (int j = 0; j < nComponents; j++) {
	if ((*mat)[i][j]) {
	  const FiniteElemSpace *feSpace = (*mat)[i][j]->getRowFeSpace();
	    
	  std::set<DegreeOfFreedom> &dirichletDofs = *((*mat)[i][j]->getApplyDBCs());

	  MSG("DIRICHLET DOFS: %d %d -> %d\n", i, j, dirichletDofs.size());
	  
	  for (std::set<DegreeOfFreedom>::iterator it = dirichletDofs.begin();
	       it != dirichletDofs.end(); ++it) {
	    if (isCoarseSpace(feSpace, *it)) {
	      if ((*coarseSpaceMap)[feSpace].isRankDof(*it)) {
		int globalDof = (*coarseSpaceMap)[feSpace][*it].global;
		dofsCoarse.push_back(coarseSpaceMap->getMatIndex(i, globalDof));
	      }
	    } else {
	      if ((*interiorMap)[feSpace].isRankDof(*it)) {
		int globalDof = (*interiorMap)[feSpace][*it].global;
		dofsInterior.push_back(interiorMap->getMatIndex(i, globalDof));
	      }
952
	    }
953
	  }
Thomas Witkowski's avatar
Thomas Witkowski committed
954
955
	} else {
	  MSG("NO MAT DIAG in %d\n", i);
956
957
958
959
960
961
962
963
964
965
	}
      }
    }

    MatZeroRows(matIntInt, dofsInterior.size(), &(dofsInterior[0]), 1.0, 
		PETSC_NULL, PETSC_NULL);

    if (coarseSpaceMap != NULL)
      MatZeroRows(matCoarseCoarse, dofsCoarse.size(), &(dofsCoarse[0]), 1.0, 
		  PETSC_NULL, PETSC_NULL);
Thomas Witkowski's avatar
Thomas Witkowski committed
966
#endif
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
  }


  void PetscSolverGlobalMatrix::removeDirichletBcDofs(SystemVector *vec)
  {
    FUNCNAME("PetscSolverGlobalMatrix::removeDirichletBcDofs()");

    int nComponents = vec->getSize();
    for (int i = 0; i < nComponents; i++) {
      const FiniteElemSpace *feSpace = vec->getDOFVector(i)->getFeSpace();

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

      for (map<DegreeOfFreedom, double>::iterator it = dirichletValues.begin();
	   it != dirichletValues.end(); ++it) {
	if (isCoarseSpace(feSpace, it->first)) {
984
985
986
	  if ((*coarseSpaceMap)[feSpace].isRankDof(it->first)) {
	    int globalDof = (*coarseSpaceMap)[feSpace][it->first].global;
	    VecSetValue(rhsCoarseSpace, coarseSpaceMap->getMatIndex(i, globalDof), 
987
			it->second, INSERT_VALUES);
988
	  }
989
990
	} else {
	  if ((*interiorMap)[feSpace].isRankDof(it->first)) {
991
992
	    int globalDof = (*interiorMap)[feSpace][it->first].global;
	    VecSetValue(rhsInterior, interiorMap->getMatIndex(i, globalDof),
993
994
995
996
997
			it->second, INSERT_VALUES);	    
	  }
	}
      }
    }
998
999
1000
1001
1002
1003
1004
1005

    VecAssemblyBegin(rhsInterior);
    VecAssemblyEnd(rhsInterior);

    if (coarseSpaceMap) {
      VecAssemblyBegin(rhsCoarseSpace);
      VecAssemblyEnd(rhsCoarseSpace);
    }
1006
1007
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
1008
}