Actual source code: ex214.c


  2: static char help[] = "Tests MatMatSolve() and MatMatTransposeSolve() for computing inv(A) with MUMPS.\n\
  3: Example: mpiexec -n <np> ./ex214 -displ \n\n";

  5: #include <petscmat.h>

  7: int main(int argc,char **args)
  8: {
  9:   PetscMPIInt    size,rank;
 10: #if defined(PETSC_HAVE_MUMPS)
 11:   Mat            A,RHS,C,F,X,AX,spRHST;
 12:   PetscInt       m,n,nrhs,M,N,i,Istart,Iend,Ii,j,J,test;
 13:   PetscScalar    v;
 14:   PetscReal      norm,tol=PETSC_SQRT_MACHINE_EPSILON;
 15:   PetscRandom    rand;
 16:   PetscBool      displ=PETSC_FALSE;
 17:   char           solver[256];
 18: #endif

 20:   PetscInitialize(&argc,&args,(char*)0,help);
 21:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 22:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

 24: #if !defined(PETSC_HAVE_MUMPS)
 25:   if (rank == 0) PetscPrintf(PETSC_COMM_SELF,"This example requires MUMPS, exit...\n");
 26:   PetscFinalize();
 27:   return 0;
 28: #else

 30:   PetscOptionsGetBool(NULL,NULL,"-displ",&displ,NULL);

 32:   /* Create matrix A */
 33:   m = 4; n = 4;
 34:   PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
 35:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);

 37:   MatCreate(PETSC_COMM_WORLD,&A);
 38:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);
 39:   MatSetFromOptions(A);
 40:   MatMPIAIJSetPreallocation(A,5,NULL,5,NULL);
 41:   MatSeqAIJSetPreallocation(A,5,NULL);

 43:   MatGetOwnershipRange(A,&Istart,&Iend);
 44:   for (Ii=Istart; Ii<Iend; Ii++) {
 45:     v = -1.0; i = Ii/n; j = Ii - i*n;
 46:     if (i>0)   {J = Ii - n; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);}
 47:     if (i<m-1) {J = Ii + n; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);}
 48:     if (j>0)   {J = Ii - 1; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);}
 49:     if (j<n-1) {J = Ii + 1; MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);}
 50:     v = 4.0; MatSetValues(A,1,&Ii,1,&Ii,&v,ADD_VALUES);
 51:   }
 52:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 53:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

 55:   MatGetLocalSize(A,&m,&n);
 56:   MatGetSize(A,&M,&N);

 59:   /* Create dense matrix C and X; C holds true solution with identical columns */
 60:   nrhs = N;
 61:   PetscOptionsGetInt(NULL,NULL,"-nrhs",&nrhs,NULL);
 62:   MatCreate(PETSC_COMM_WORLD,&C);
 63:   MatSetSizes(C,m,PETSC_DECIDE,PETSC_DECIDE,nrhs);
 64:   MatSetType(C,MATDENSE);
 65:   MatSetFromOptions(C);
 66:   MatSetUp(C);

 68:   PetscRandomCreate(PETSC_COMM_WORLD,&rand);
 69:   PetscRandomSetFromOptions(rand);
 70:   MatSetRandom(C,rand);
 71:   MatDuplicate(C,MAT_DO_NOT_COPY_VALUES,&X);

 73:   PetscStrcpy(solver,MATSOLVERMUMPS);
 74:   if (rank == 0 && displ) PetscPrintf(PETSC_COMM_SELF,"Solving with %s: nrhs %" PetscInt_FMT ", size mat %" PetscInt_FMT " x %" PetscInt_FMT "\n",solver,nrhs,M,N);

 76:   for (test=0; test<2; test++) {
 77:     if (test == 0) {
 78:       /* Test LU Factorization */
 79:       PetscPrintf(PETSC_COMM_WORLD,"test LU factorization\n");
 80:       MatGetFactor(A,solver,MAT_FACTOR_LU,&F);
 81:       MatLUFactorSymbolic(F,A,NULL,NULL,NULL);
 82:       MatLUFactorNumeric(F,A,NULL);
 83:     } else {
 84:       /* Test Cholesky Factorization */
 85:       PetscBool flg;
 86:       MatIsSymmetric(A,0.0,&flg);

 89:       PetscPrintf(PETSC_COMM_WORLD,"test Cholesky factorization\n");
 90:       MatGetFactor(A,solver,MAT_FACTOR_CHOLESKY,&F);
 91:       MatCholeskyFactorSymbolic(F,A,NULL,NULL);
 92:       MatCholeskyFactorNumeric(F,A,NULL);
 93:     }

 95:     /* (1) Test MatMatSolve(): dense RHS = A*C, C: true solutions */
 96:     /* ---------------------------------------------------------- */
 97:     MatMatMult(A,C,MAT_INITIAL_MATRIX,2.0,&RHS);
 98:     MatMatSolve(F,RHS,X);
 99:     /* Check the error */
100:     MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);
101:     MatNorm(X,NORM_FROBENIUS,&norm);
102:     if (norm > tol) {
103:       PetscPrintf(PETSC_COMM_SELF,"(1) MatMatSolve: Norm of error %g\n",norm);
104:     }

106:     /* Test X=RHS */
107:     MatMatSolve(F,RHS,RHS);
108:     /* Check the error */
109:     MatAXPY(RHS,-1.0,C,SAME_NONZERO_PATTERN);
110:     MatNorm(RHS,NORM_FROBENIUS,&norm);
111:     if (norm > tol) {
112:       PetscPrintf(PETSC_COMM_SELF,"(1.1) MatMatSolve: Norm of error %g\n",norm);
113:     }

115:     /* (2) Test MatMatSolve() for inv(A) with dense RHS:
116:      RHS = [e[0],...,e[nrhs-1]], dense X holds first nrhs columns of inv(A) */
117:     /* -------------------------------------------------------------------- */
118:     MatZeroEntries(RHS);
119:     for (i=0; i<nrhs; i++) {
120:       v = 1.0;
121:       MatSetValues(RHS,1,&i,1,&i,&v,INSERT_VALUES);
122:     }
123:     MatAssemblyBegin(RHS,MAT_FINAL_ASSEMBLY);
124:     MatAssemblyEnd(RHS,MAT_FINAL_ASSEMBLY);

126:     MatMatSolve(F,RHS,X);
127:     if (displ) {
128:       if (rank == 0) PetscPrintf(PETSC_COMM_SELF," \n(2) first %" PetscInt_FMT " columns of inv(A) with dense RHS:\n",nrhs);
129:       MatView(X,PETSC_VIEWER_STDOUT_WORLD);
130:     }

132:     /* Check the residual */
133:     MatMatMult(A,X,MAT_INITIAL_MATRIX,2.0,&AX);
134:     MatAXPY(AX,-1.0,RHS,SAME_NONZERO_PATTERN);
135:     MatNorm(AX,NORM_INFINITY,&norm);
136:     if (norm > tol) {
137:       PetscPrintf(PETSC_COMM_SELF,"(2) MatMatSolve: Norm of residual %g\n",norm);
138:     }
139:     MatZeroEntries(X);

141:     /* (3) Test MatMatTransposeSolve() for inv(A) with sparse RHS stored in the host:
142:      spRHST = [e[0],...,e[nrhs-1]]^T, dense X holds first nrhs columns of inv(A) */
143:     /* --------------------------------------------------------------------------- */
144:     /* Create spRHST: PETSc does not support compressed column format which is required by MUMPS for sparse RHS matrix,
145:      thus user must create spRHST=spRHS^T and call MatMatTransposeSolve() */
146:     MatCreate(PETSC_COMM_WORLD,&spRHST);
147:     if (rank == 0) {
148:       /* MUMPS requirs RHS be centralized on the host! */
149:       MatSetSizes(spRHST,nrhs,M,PETSC_DECIDE,PETSC_DECIDE);
150:     } else {
151:       MatSetSizes(spRHST,0,0,PETSC_DECIDE,PETSC_DECIDE);
152:     }
153:     MatSetType(spRHST,MATAIJ);
154:     MatSetFromOptions(spRHST);
155:     MatSetUp(spRHST);
156:     if (rank == 0) {
157:       v = 1.0;
158:       for (i=0; i<nrhs; i++) {
159:         MatSetValues(spRHST,1,&i,1,&i,&v,INSERT_VALUES);
160:       }
161:     }
162:     MatAssemblyBegin(spRHST,MAT_FINAL_ASSEMBLY);
163:     MatAssemblyEnd(spRHST,MAT_FINAL_ASSEMBLY);

165:     MatMatTransposeSolve(F,spRHST,X);

167:     if (displ) {
168:       if (rank == 0) PetscPrintf(PETSC_COMM_SELF," \n(3) first %" PetscInt_FMT " columns of inv(A) with sparse RHS:\n",nrhs);
169:       MatView(X,PETSC_VIEWER_STDOUT_WORLD);
170:     }

172:     /* Check the residual: R = A*X - RHS */
173:     MatMatMult(A,X,MAT_REUSE_MATRIX,2.0,&AX);

175:     MatAXPY(AX,-1.0,RHS,SAME_NONZERO_PATTERN);
176:     MatNorm(AX,NORM_INFINITY,&norm);
177:     if (norm > tol) {
178:       PetscPrintf(PETSC_COMM_SELF,"(3) MatMatSolve: Norm of residual %g\n",norm);
179:     }

181:     /* (4) Test MatMatSolve() for inv(A) with selected entries:
182:      input: spRHS gives selected indices; output: spRHS holds selected entries of inv(A) */
183:     /* --------------------------------------------------------------------------------- */
184:     if (nrhs == N) { /* mumps requires nrhs = n */
185:       /* Create spRHS on proc[0] */
186:       Mat spRHS = NULL;

188:       /* Create spRHS = spRHST^T. Two matrices share internal matrix data structure */
189:       MatCreateTranspose(spRHST,&spRHS);
190:       MatMumpsGetInverse(F,spRHS);
191:       MatDestroy(&spRHS);

193:       MatMumpsGetInverseTranspose(F,spRHST);
194:       if (displ) {
195:         PetscPrintf(PETSC_COMM_WORLD,"\nSelected entries of inv(A^T):\n");
196:         MatView(spRHST,PETSC_VIEWER_STDOUT_WORLD);
197:       }
198:       MatDestroy(&spRHS);
199:     }
200:     MatDestroy(&AX);
201:     MatDestroy(&F);
202:     MatDestroy(&RHS);
203:     MatDestroy(&spRHST);
204:   }

206:   /* Free data structures */
207:   MatDestroy(&A);
208:   MatDestroy(&C);
209:   MatDestroy(&X);
210:   PetscRandomDestroy(&rand);
211:   PetscFinalize();
212:   return 0;
213: #endif
214: }

216: /*TEST

218:    test:
219:      requires: mumps double !complex

221:    test:
222:      suffix: 2
223:      requires: mumps double !complex
224:      nsize: 2

226: TEST*/