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

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

24
25
26
27
    petscData.init(interiorMap, coarseSpaceMap, 
		   subdomainLevel, mpiCommLocal, mpiCommGlobal,
		   meshDistributor);
    petscData.create();
28

Thomas Witkowski's avatar
Thomas Witkowski committed
29
    if (coarseSpaceMap.size()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
30
      updateSubdomainData();
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
31
      fillPetscMatrixWithCoarseSpace(seqMat);
32
33
34
      return;
    }

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

41
42
    // === Create PETSc vector (solution and a temporary vector). ===

43
44
    int nRankRows = interiorMap->getRankDofs();
    int nOverallRows = interiorMap->getOverallDofs();
45

46
    VecCreateMPI(mpiCommGlobal, nRankRows, nOverallRows, &petscSolVec);
47
48


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

#if (DEBUG != 0)
    int a, b;
55
    MatGetOwnershipRange(petscData.getInteriorMat(), &a, &b);
56
57
    TEST_EXIT(a == interiorMap->getStartDofs())("Wrong matrix ownership range!\n");
    TEST_EXIT(b == interiorMap->getStartDofs() + nRankRows)
Thomas Witkowski's avatar
Thomas Witkowski committed
58
59
60
61
62
63
      ("Wrong matrix ownership range!\n");
#endif


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

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
64
    int nComponents = seqMat->getNumRows();
Thomas Witkowski's avatar
Thomas Witkowski committed
65
66
    for (int i = 0; i < nComponents; i++)
      for (int j = 0; j < nComponents; j++)
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
67
68
	if ((*seqMat)[i][j])
	  setDofMatrix((*seqMat)[i][j], i, j);
Thomas Witkowski's avatar
Thomas Witkowski committed
69
70
71
72
73

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

74
    petscData.assembly();
75
76
77
    
    if (printMatInfo) {
      MatInfo matInfo;
78
      MatGetInfo(petscData.getInteriorMat(), MAT_GLOBAL_SUM, &matInfo);
79
80
81
82
83
84
85
      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));
    }
86
87
88

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

89
    //    removeDirichletBcDofs(mat);
90
91
    

92
    // === Init PETSc solver. ===
93

94
95
    KSPCreate(mpiCommGlobal, &kspInterior);
    KSPGetPC(kspInterior, &pcInterior);
96
97
98
99
    KSPSetOperators(kspInterior, 
		    petscData.getInteriorMat(), 
		    petscData.getInteriorMat(), 
		    SAME_NONZERO_PATTERN); 
100
101
102
103
    KSPSetTolerances(kspInterior, 0.0, 1e-8, PETSC_DEFAULT, PETSC_DEFAULT);
    KSPSetType(kspInterior, KSPBCGS);
    KSPSetOptionsPrefix(kspInterior, kspPrefix.c_str());
    KSPSetFromOptions(kspInterior);
104

105
    initPreconditioner(pcInterior);
106

107
108
    // Do not delete the solution vector, use it for the initial guess.
    if (!zeroStartVector)
109
      KSPSetInitialGuessNonzero(kspInterior, PETSC_TRUE);
110

Thomas Witkowski's avatar
Thomas Witkowski committed
111
112
113
#if (DEBUG != 0)
    MSG("Fill petsc matrix 3 needed %.5f seconds\n", MPI::Wtime() - wtime);
#endif
114
115
116
  }


117
118
119
120
121
122
123
124
  void PetscSolverGlobalMatrix::fillPetscMatrix(DOFMatrix *mat)
  {
    Matrix<DOFMatrix*> m(1, 1);
    m[0][0] = mat;
    fillPetscMatrix(&m);
  }


Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
125
  void PetscSolverGlobalMatrix::fillPetscMatrixWithCoarseSpace(Matrix<DOFMatrix*> *seqMat)
126
127
  {
    FUNCNAME("PetscSolverGlobalMatrix::fillPetscMatrixWithCoarseSpace()");
128
129

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

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
133
134
    vector<const FiniteElemSpace*> feSpaces = getFeSpaces(seqMat);

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

138

139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
    // === 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
158
    int nComponents = seqMat->getSize();
159
160
    for (int i = 0; i < nComponents; i++) {
      for (int j = 0; j < nComponents; j++) {
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
161
	if (!(*seqMat)[i][j])
162
163
	  continue;

Thomas Witkowski's avatar
Thomas Witkowski committed
164
165
166
	ParallelDofMapping *rowCoarseSpace = coarseSpaceMap[i];
	ParallelDofMapping *colCoarseSpace = coarseSpaceMap[j];

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
167
168
	traits::col<Matrix>::type col((*seqMat)[i][j]->getBaseMatrix());
	traits::const_value<Matrix>::type value((*seqMat)[i][j]->getBaseMatrix());
169
170
	
	// Traverse all rows.
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
171
172
	for (cursor_type cursor = begin<row>((*seqMat)[i][j]->getBaseMatrix()), 
	       cend = end<row>((*seqMat)[i][j]->getBaseMatrix()); cursor != cend; ++cursor) {
173

Thomas Witkowski's avatar
Thomas Witkowski committed
174
	  bool isRowCoarse = isCoarseSpace(i, feSpaces[i], *cursor);
175
176
177
178
179
180
181
182
183
184
  
	  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
185
	    bool isColCoarse = isCoarseSpace(j, feSpaces[j], col(*icursor));
186

Thomas Witkowski's avatar
Thomas Witkowski committed
187
188
	    if (isColCoarse) {
	      if (isRowCoarse) {
189
190
191
192
193
194
195
		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
196
	      if (isRowCoarse) {
197
198
199
200
201
202
203
204
205
206
207
208
		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
209
210
	  if (isRowCoarse) {
	    int rowIndex = rowCoarseSpace->getMatIndex(i, *cursor);
211
	    for (unsigned int k = 0; k < cols.size(); k++)
Thomas Witkowski's avatar
Thomas Witkowski committed
212
	      cols[k] = colCoarseSpace->getMatIndex(j, cols[k]);
213

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
214
	    // matCoarseCoarse
215
	    MatSetValues(petscData.getCoarseMatComp(i),
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
216
			 1, &rowIndex, cols.size(),
217
218
219
220
			 &(cols[0]), &(values[0]), ADD_VALUES);

	    if (colsOther.size()) {
	      if (subdomainLevel == 0) {
221
		for (unsigned int k = 0; k < colsOther.size(); k++)
222
223
224
225
226
227
		  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;
	      }
228

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
229
	      // matCoarseInt
230
	      MatSetValues(petscData.getCoarseIntMatComp(i), 
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
231
			   1, &rowIndex, colsOther.size(),
232
233
234
235
236
237
238
239
240
241
242
243
244
245
 			   &(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]);
	    }
	    
246
  	    MatSetValues(petscData.getInteriorMat(), 1, &localRowIndex, cols.size(),
247
248
249
250
251
252
253
254
255
  			 &(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++)
Thomas Witkowski's avatar
Thomas Witkowski committed
256
		colsOther[k] = colCoarseSpace->getMatIndex(j, colsOther[k]);
257

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
258
	      // matIntCoarse
259
  	      MatSetValues(petscData.getIntCoarseMatComp(i), 
Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
260
			   1, &globalRowIndex, colsOther.size(),
261
262
263
264
265
266
267
  			   &(colsOther[0]), &(valuesOther[0]), ADD_VALUES);
	    }
	  }
	} 
      }
    }

268
    petscData.assembly();
269

270
271
272
    // === Create solver for the non primal (thus local) variables. ===

    KSPCreate(mpiCommLocal, &kspInterior);
273
274
275
276
    KSPSetOperators(kspInterior, 
		    petscData.getInteriorMat(),
		    petscData.getInteriorMat(),
		    SAME_NONZERO_PATTERN);
277
278
279
280
281
282
283
284
285
286
287
288
    KSPSetOptionsPrefix(kspInterior, "interior_");
    KSPSetType(kspInterior, KSPPREONLY);
    KSPGetPC(kspInterior, &pcInterior);
    PCSetType(pcInterior, PCLU);
    if (subdomainLevel == 0)
      PCFactorSetMatSolverPackage(pcInterior, MATSOLVERUMFPACK);
    else
      PCFactorSetMatSolverPackage(pcInterior, MATSOLVERMUMPS);
    KSPSetFromOptions(kspInterior);  
  }


289
290
291
  void PetscSolverGlobalMatrix::fillPetscRhs(SystemVector *vec)
  {
    FUNCNAME("PetscSolverGlobalMatrix::fillPetscRhs()");
Thomas Witkowski's avatar
Thomas Witkowski committed
292

Thomas Witkowski's avatar
Thomas Witkowski committed
293
294
295
296
    VecCreateMPI(mpiCommGlobal, 
		 interiorMap->getRankDofs(), 
		 nGlobalOverallInterior,
		 &rhsInterior);
297

Thomas Witkowski's avatar
Thomas Witkowski committed
298
    if (coarseSpaceMap.size()) 
Thomas Witkowski's avatar
Thomas Witkowski committed
299
      VecCreateMPI(mpiCommGlobal, 
Thomas Witkowski's avatar
Thomas Witkowski committed
300
301
		   coarseSpaceMap[0]->getRankDofs(), 
		   coarseSpaceMap[0]->getOverallDofs(),
Thomas Witkowski's avatar
Thomas Witkowski committed
302
303
		   &rhsCoarseSpace);
    
304
    TEST_EXIT_DBG(vec)("No DOF vector defined!\n");
305
    TEST_EXIT_DBG(interiorMap)("No parallel DOF map defined!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
306
    
307

Thomas Witkowski's avatar
Thomas Witkowski committed
308
    // === Transfer values from DOF vector to the PETSc vector. === 
Thomas Witkowski's avatar
Thomas Witkowski committed
309
    if (coarseSpaceMap.size()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
310
311
      for (int i = 0; i < vec->getSize(); i++)
	setDofVector(rhsInterior, rhsCoarseSpace, vec->getDOFVector(i), i);
Thomas Witkowski's avatar
Thomas Witkowski committed
312
313
314
315
    } else {
      for (int i = 0; i < vec->getSize(); i++)
	setDofVector(rhsInterior, vec->getDOFVector(i), i);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
316

317
318
    VecAssemblyBegin(rhsInterior);
    VecAssemblyEnd(rhsInterior);
319

Thomas Witkowski's avatar
Thomas Witkowski committed
320
    if (coarseSpaceMap.size()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
321
322
323
324
      VecAssemblyBegin(rhsCoarseSpace);
      VecAssemblyEnd(rhsCoarseSpace);
    }

325
326
327

    // === Remove Dirichlet BC DOFs. ===
    //    removeDirichletBcDofs(vec);
Thomas Witkowski's avatar
Thomas Witkowski committed
328
 
Thomas Witkowski's avatar
Thomas Witkowski committed
329
330
331
  }


332
333
  void PetscSolverGlobalMatrix::solvePetscMatrix(SystemVector &vec, 
						 AdaptInfo *adaptInfo)
Thomas Witkowski's avatar
Thomas Witkowski committed
334
335
336
337
338
339
340
341
342
343
  {
    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++)
344
	setDofVector(petscSolVec, vec.getDOFVector(i), i, true);
Thomas Witkowski's avatar
Thomas Witkowski committed
345
346
347
348
349
      
      VecAssemblyBegin(petscSolVec);
      VecAssemblyEnd(petscSolVec);
    }

350

Thomas Witkowski's avatar
Thomas Witkowski committed
351
352
    MatNullSpace matNullspace;
    Vec nullspaceBasis;
353
354
355
    if (nullspace.size() > 0 || 
	hasConstantNullspace ||
	constNullspaceComponent.size() > 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
356
      TEST_EXIT_DBG(nullspace.size() <= 1)("Not yet implemented!\n");
357

358
359
360
361
362
363
364
365
366
367
      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
368
369
370
371
372
373
374
      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);
375
376

	VecNormalize(nullspaceBasis, PETSC_NULL);
Thomas Witkowski's avatar
Thomas Witkowski committed
377
378
379
380
	
	MatNullSpaceCreate(mpiCommGlobal, (hasConstantNullspace ? PETSC_TRUE : PETSC_FALSE), 
			   1, &nullspaceBasis, &matNullspace);

381
	MatMult(petscData.getInteriorMat(), nullspaceBasis, petscSolVec);
Thomas Witkowski's avatar
Thomas Witkowski committed
382
383
384
385
386
387
	PetscReal n;
	VecNorm(petscSolVec, NORM_2, &n);
	MSG("NORM IS: %e\n", n);
      } else {
	MatNullSpaceCreate(mpiCommGlobal, PETSC_TRUE, 0, PETSC_NULL, &matNullspace);
      }
388

389
      MatSetNullSpace(petscData.getInteriorMat(), matNullspace);
Thomas Witkowski's avatar
Thomas Witkowski committed
390
      KSPSetNullSpace(kspInterior, matNullspace);
391

Thomas Witkowski's avatar
Thomas Witkowski committed
392
      // === Remove null space, if requested. ===
393

Thomas Witkowski's avatar
Thomas Witkowski committed
394
      if (removeRhsNullspace) {
Thomas Witkowski's avatar
Thomas Witkowski committed
395
	TEST_EXIT_DBG(coarseSpaceMap.empty())("Not supported!\n");
396
397

	MSG("Remove nullspace from rhs vector.\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
398
399
400
401
402
403
	
	MatNullSpaceRemove(matNullspace, rhsInterior, PETSC_NULL);
      }
    } else {
      TEST_EXIT(removeRhsNullspace == false)
	("No nullspace provided that should be removed from rhs!\n");
404
405
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
406
    // PETSc.
407
    solve(rhsInterior, petscSolVec);
Thomas Witkowski's avatar
Thomas Witkowski committed
408

409
410

    if (nullspace.size() > 0) {
Thomas Witkowski's avatar
Thomas Witkowski committed
411
      MatNullSpaceDestroy(&matNullspace);
412
413
414
415
      VecDestroy(&nullspaceBasis);
    }


Thomas Witkowski's avatar
Thomas Witkowski committed
416
417
    // === Transfere values from PETSc's solution vectors to the DOF vectors. ===
    PetscScalar *vecPointer;
418
    VecGetArray(petscSolVec, &vecPointer);    
Thomas Witkowski's avatar
Thomas Witkowski committed
419

420
    int c = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
421
    for (int i = 0; i < nComponents; i++) {
422
      DOFVector<double> &dv = *(vec.getDOFVector(i));
423

424
      DofMap& d = (*interiorMap)[dv.getFeSpace()].getMap();
425
426
427
      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
428
429
430
    }

    VecRestoreArray(petscSolVec, &vecPointer);
Thomas Witkowski's avatar
Thomas Witkowski committed
431
432
433
434

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


438
439
440
441
  void PetscSolverGlobalMatrix::solveGlobal(Vec &rhs, Vec &sol)
  {
    FUNCNAME("PetscSolverGlobalMatrix::solveGlobal()");

442
443
444
    double wtime = MPI::Wtime();
    double t0 = 0.0, t1 = 0.0;

445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
    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);

464
465
466
    t0 = MPI::Wtime() - wtime;

    wtime = MPI::Wtime();
467
    KSPSolve(kspInterior, tmp, tmp);
468
469
470
    t1 = MPI::Wtime() - wtime;

    wtime = MPI::Wtime();
471
472
473
474
475
476
477
478
479
480
481

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

484
    //    MSG("TIMEING: %.5f %.5f\n", t0, t1);
485
486
487
  }


488
489
490
491
  void PetscSolverGlobalMatrix::destroyMatrixData()
  {
    FUNCNAME("PetscSolverGlobalMatrix::destroyMatrixData()");

492
493
    exitPreconditioner(pcInterior);

494
    petscData.destroy();
Thomas Witkowski's avatar
Thomas Witkowski committed
495

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
496
    KSPDestroy(&kspInterior);
497
498

    if (petscSolVec != PETSC_NULL) {
Thomas Witkowski's avatar
Thomas Witkowski committed
499
      VecDestroy(&petscSolVec);
500
501
      petscSolVec = PETSC_NULL;
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
502
503
504
  }


505
506
507
508
  void PetscSolverGlobalMatrix::destroyVectorData()
  {
    FUNCNAME("PetscSolverGlobalMatrix::destroyVectorData()");

509
    VecDestroy(&rhsInterior);
510

Thomas Witkowski's avatar
Thomas Witkowski committed
511
    if (coarseSpaceMap.size())
512
      VecDestroy(&rhsCoarseSpace);
513
514
515
  }


516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
  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);
    }
  }


552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
  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
567
  void PetscSolverGlobalMatrix::setDofMatrix(DOFMatrix* seqMat,
568
					     int nRowMat, int nColMat)
Thomas Witkowski's avatar
Thomas Witkowski committed
569
570
571
  {
    FUNCNAME("PetscSolverGlobalMatrix::setDofMatrix()");

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
572
    TEST_EXIT(seqMat)("No DOFMatrix!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
573
574

    using mtl::tag::row; using mtl::tag::nz; using mtl::begin; using mtl::end;
575
    namespace traits = mtl::traits;
Thomas Witkowski's avatar
Thomas Witkowski committed
576
577
    typedef DOFMatrix::base_matrix_type Matrix;

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
578
579
    traits::col<Matrix>::type col(seqMat->getBaseMatrix());
    traits::const_value<Matrix>::type value(seqMat->getBaseMatrix());
Thomas Witkowski's avatar
Thomas Witkowski committed
580
581
582
583
584
585
586
587
588
589
590

    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;

591
592
593
    // Get periodic mapping object
    PeriodicMap &perMap = meshDistributor->getPeriodicMap();

Thomas Witkowski's avatar
Blub  
Thomas Witkowski committed
594
595
    const FiniteElemSpace *rowFe = seqMat->getRowFeSpace();
    const FiniteElemSpace *colFe = seqMat->getColFeSpace();
Thomas Witkowski's avatar
Thomas Witkowski committed
596
597
598
599
    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
600
601
    // === to the PETSc matrix.                                               ===

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
610
611
612
      if (!periodicRow) {
	// === Row DOF index is not periodic. ===

613
	// Get PETSc's mat row index.
614
	int rowIndex = interiorMap->getMatIndex(nRowMat, globalRowDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
615
616
617
618
619
620
621
622

	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
623
	  int globalColDof = colMap[col(*icursor)].global;
Thomas Witkowski's avatar
Thomas Witkowski committed
624
	  // Test if the current col dof is a periodic dof.
625
	  bool periodicCol = perMap.isPeriodic(colFe, globalColDof);
626
	  // Get PETSc's mat col index.
627
	  int colIndex = interiorMap->getMatIndex(nColMat, globalColDof);
Thomas Witkowski's avatar
Thomas Witkowski committed
628
629
630
631
632
633
634

	  // 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
635
636
 	    cols.push_back(colIndex);
 	    values.push_back(value(*icursor));
Thomas Witkowski's avatar
Thomas Witkowski committed
637
638
639
640
641
	  } else {
	    // === Row index is not periodic, but column index is. ===

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

Thomas Witkowski's avatar
Thomas Witkowski committed
645
646
647
648
649
650
651
652
653
	    // 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;
654
	    perMap.mapDof(colFe, globalColDof, perAsc, newCols);
Thomas Witkowski's avatar
Thomas Witkowski committed
655
	    for (unsigned int i = 0; i < newCols.size(); i++) {
656
	      cols.push_back(interiorMap->getMatIndex(nColMat, newCols[i]));
Thomas Witkowski's avatar
Thomas Witkowski committed
657
658
659
660
661
	      values.push_back(scaledValue);	      
	    }
	  }
	}

662
  	MatSetValues(petscData.getInteriorMat(), 1, &rowIndex, cols.size(), 
Thomas Witkowski's avatar
Thomas Witkowski committed
663
  		     &(cols[0]), &(values[0]), ADD_VALUES);	
Thomas Witkowski's avatar
Thomas Witkowski committed
664
665
666
667
668
669
670
671
672
673
674
675
676
677
      } 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.
678
	  int globalColDof = (*interiorMap)[colFe][col(*icursor)].global;
Thomas Witkowski's avatar
Thomas Witkowski committed
679
680
681
682
683
684
685
686
687

	  // 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;
688
689
690
691
	  perMap.fillAssociations(colFe, globalColDof, 
				  meshDistributor->getElementObjectDb(), perAsc);
	  perMap.fillAssociations(rowFe, globalRowDof, 
				  meshDistributor->getElementObjectDb(), perAsc);
Thomas Witkowski's avatar
Thomas Witkowski committed
692
693
694
695
696
697
698
699
700
701

	  // 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
702
703
	  perMap.mapDof(rowFe, colFe, make_pair(globalRowDof, globalColDof),
			perAsc, entry);
Thomas Witkowski's avatar
Thomas Witkowski committed
704
705
706

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

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

711
712
	    colsMap[rowIdx].push_back(colIdx);
	    valsMap[rowIdx].push_back(scaledValue);
Thomas Witkowski's avatar
Thomas Witkowski committed
713
714
715
716
717
718
719
720
721
722
723
724
	  }
	}


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

726
	  MatSetValues(petscData.getInteriorMat(), 1, &rowIndex, rowIt->second.size(),
Thomas Witkowski's avatar
Thomas Witkowski committed
727
728
729
730
731
732
733
		       &(rowIt->second[0]), &(valsMap[rowIt->first][0]), ADD_VALUES);
	}
      }
    }
  }


Thomas Witkowski's avatar
Thomas Witkowski committed
734
735
  void PetscSolverGlobalMatrix::setDofVector(Vec vecInterior, 
					     Vec vecCoarse,
736
					     DOFVector<double>* vec, 
737
					     int nRowVec, 
738
					     bool rankOnly)
Thomas Witkowski's avatar
Thomas Witkowski committed
739
740
741
  {
    FUNCNAME("PetscSolverGlobalMatrix::setDofVector()");

742
    const FiniteElemSpace *feSpace = vec->getFeSpace();
743
    PeriodicMap &perMap = meshDistributor->getPeriodicMap();
Thomas Witkowski's avatar
Thomas Witkowski committed
744
745
746
747
    
    ParallelDofMapping *rowCoarseSpace = NULL;
    if (coarseSpaceMap.size())
      rowCoarseSpace = coarseSpaceMap[nRowVec];
748

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

753
      if (rankOnly && !(*interiorMap)[feSpace].isRankDof(dofIt.getDOFIndex()))
Thomas Witkowski's avatar
Thomas Witkowski committed
754
755
	continue;

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

Thomas Witkowski's avatar
Thomas Witkowski committed
759
	int index = rowCoarseSpace->getMatIndex(nRowVec, dofIt.getDOFIndex());
760
761
762
763
764
765
766
767
768
769
770
771
772
773
	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;
774

775
	if (perMap.isPeriodic(feSpace, globalRowDof)) {
Thomas Witkowski's avatar
Thomas Witkowski committed
776
777
778
	  std::set<int>& perAsc = perMap.getAssociations(feSpace, globalRowDof);
	  double value = *dofIt / (perAsc.size() + 1.0);
	  VecSetValue(vecInterior, index, value, ADD_VALUES);
779

Thomas Witkowski's avatar
Thomas Witkowski committed
780
781
782
783
	  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);
784

Thomas Witkowski's avatar
Thomas Witkowski committed
785
	    VecSetValue(vecInterior, mappedIndex, value, ADD_VALUES);
786
787
788
789
	  }	  
	} else {	  
	  // The DOF index is not periodic.
	  
Thomas Witkowski's avatar
Thomas Witkowski committed
790
791
	  VecSetValue(vecInterior, index, *dofIt, ADD_VALUES);
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
792
793
794
795
796
      }
    }
  }


797
798
  void PetscSolverGlobalMatrix::removeDirichletBcDofs(Matrix<DOFMatrix*> *mat)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
799
    FUNCNAME("PetscSolverGlobalMatrix::removeDirichletBcDofs()");    
800
801
802
803
804
805
806
807
  }


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

Thomas Witkowski's avatar
Thomas Witkowski committed
808
}