Actual source code: ex245.c


  2: static char help[] = "Tests LU, Cholesky factorization and MatMatSolve() for a ScaLAPACK dense matrix.\n\n";

  4: #include <petscmat.h>

  6: int main(int argc,char **argv)
  7: {
  8:   Mat            A,F,B,X,C,Aher,G;
  9:   Vec            b,x,c,d,e;
 11:   PetscInt       m=5,n,p,i,j,nrows,ncols;
 12:   PetscScalar    *v,*barray,rval;
 13:   PetscReal      norm,tol=1.e5*PETSC_MACHINE_EPSILON;
 14:   PetscMPIInt    size,rank;
 15:   PetscRandom    rand;
 16:   const PetscInt *rows,*cols;
 17:   IS             isrows,iscols;
 18:   PetscBool      mats_view=PETSC_FALSE;

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

 24:   PetscRandomCreate(PETSC_COMM_WORLD,&rand);
 25:   PetscRandomSetFromOptions(rand);

 27:   /* Get local dimensions of matrices */
 28:   PetscOptionsGetInt(NULL,NULL,"-m",&m,NULL);
 29:   n    = m;
 30:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
 31:   p    = m/2;
 32:   PetscOptionsGetInt(NULL,NULL,"-p",&p,NULL);
 33:   PetscOptionsHasName(NULL,NULL,"-mats_view",&mats_view);

 35:   /* Create matrix A */
 36:   PetscPrintf(PETSC_COMM_WORLD," Create ScaLAPACK matrix A\n");
 37:   MatCreate(PETSC_COMM_WORLD,&A);
 38:   MatSetSizes(A,m,n,PETSC_DECIDE,PETSC_DECIDE);
 39:   MatSetType(A,MATSCALAPACK);
 40:   MatSetFromOptions(A);
 41:   MatSetUp(A);
 42:   /* Set local matrix entries */
 43:   MatGetOwnershipIS(A,&isrows,&iscols);
 44:   ISGetLocalSize(isrows,&nrows);
 45:   ISGetIndices(isrows,&rows);
 46:   ISGetLocalSize(iscols,&ncols);
 47:   ISGetIndices(iscols,&cols);
 48:   PetscMalloc1(nrows*ncols,&v);
 49:   for (i=0;i<nrows;i++) {
 50:     for (j=0;j<ncols;j++) {
 51:       PetscRandomGetValue(rand,&rval);
 52:       v[i*ncols+j] = rval;
 53:     }
 54:   }
 55:   MatSetValues(A,nrows,rows,ncols,cols,v,INSERT_VALUES);
 56:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 57:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
 58:   ISRestoreIndices(isrows,&rows);
 59:   ISRestoreIndices(iscols,&cols);
 60:   ISDestroy(&isrows);
 61:   ISDestroy(&iscols);
 62:   PetscFree(v);
 63:   if (mats_view) {
 64:     PetscPrintf(PETSC_COMM_WORLD, "A: nrows %d, m %d; ncols %d, n %d\n",nrows,m,ncols,n);
 65:     MatView(A,PETSC_VIEWER_STDOUT_WORLD);
 66:   }

 68:   /* Create rhs matrix B */
 69:   PetscPrintf(PETSC_COMM_WORLD," Create rhs matrix B\n");
 70:   MatCreate(PETSC_COMM_WORLD,&B);
 71:   MatSetSizes(B,m,p,PETSC_DECIDE,PETSC_DECIDE);
 72:   MatSetType(B,MATSCALAPACK);
 73:   MatSetFromOptions(B);
 74:   MatSetUp(B);
 75:   MatGetOwnershipIS(B,&isrows,&iscols);
 76:   ISGetLocalSize(isrows,&nrows);
 77:   ISGetIndices(isrows,&rows);
 78:   ISGetLocalSize(iscols,&ncols);
 79:   ISGetIndices(iscols,&cols);
 80:   PetscMalloc1(nrows*ncols,&v);
 81:   for (i=0;i<nrows;i++) {
 82:     for (j=0;j<ncols;j++) {
 83:       PetscRandomGetValue(rand,&rval);
 84:       v[i*ncols+j] = rval;
 85:     }
 86:   }
 87:   MatSetValues(B,nrows,rows,ncols,cols,v,INSERT_VALUES);
 88:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
 89:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
 90:   ISRestoreIndices(isrows,&rows);
 91:   ISRestoreIndices(iscols,&cols);
 92:   ISDestroy(&isrows);
 93:   ISDestroy(&iscols);
 94:   PetscFree(v);
 95:   if (mats_view) {
 96:     PetscPrintf(PETSC_COMM_WORLD, "B: nrows %d, m %d; ncols %d, p %d\n",nrows,m,ncols,p);
 97:     MatView(B,PETSC_VIEWER_STDOUT_WORLD);
 98:   }

100:   /* Create rhs vector b and solution x (same size as b) */
101:   VecCreate(PETSC_COMM_WORLD,&b);
102:   VecSetSizes(b,m,PETSC_DECIDE);
103:   VecSetFromOptions(b);
104:   VecGetArray(b,&barray);
105:   for (j=0;j<m;j++) {
106:     PetscRandomGetValue(rand,&rval);
107:     barray[j] = rval;
108:   }
109:   VecRestoreArray(b,&barray);
110:   VecAssemblyBegin(b);
111:   VecAssemblyEnd(b);
112:   if (mats_view) {
113:     PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d] b: m %d\n",rank,m);
114:     PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
115:     VecView(b,PETSC_VIEWER_STDOUT_WORLD);
116:   }
117:   VecDuplicate(b,&x);

119:   /* Create matrix X - same size as B */
120:   PetscPrintf(PETSC_COMM_WORLD," Create solution matrix X\n");
121:   MatDuplicate(B,MAT_DO_NOT_COPY_VALUES,&X);

123:   /* Cholesky factorization */
124:   /*------------------------*/
125:   PetscPrintf(PETSC_COMM_WORLD," Create ScaLAPACK matrix Aher\n");
126:   MatHermitianTranspose(A,MAT_INITIAL_MATRIX,&Aher);
127:   MatAXPY(Aher,1.0,A,SAME_NONZERO_PATTERN); /* Aher = A + A^T */
128:   MatShift(Aher,100.0);  /* add 100.0 to diagonals of Aher to make it spd */
129:   if (mats_view) {
130:     PetscPrintf(PETSC_COMM_WORLD, "Aher:\n");
131:     MatView(Aher,PETSC_VIEWER_STDOUT_WORLD);
132:   }

134:   /* Cholesky factorization */
135:   /*------------------------*/
136:   PetscPrintf(PETSC_COMM_WORLD," Test Cholesky Solver \n");
137:   /* In-place Cholesky */
138:   /* Create matrix factor G, with a copy of Aher */
139:   MatDuplicate(Aher,MAT_COPY_VALUES,&G);

141:   /* G = L * L^T */
142:   MatCholeskyFactor(G,0,0);
143:   if (mats_view) {
144:     PetscPrintf(PETSC_COMM_WORLD, "Cholesky Factor G:\n");
145:     MatView(G,PETSC_VIEWER_STDOUT_WORLD);
146:   }

148:   /* Solve L * L^T x = b and L * L^T * X = B */
149:   MatSolve(G,b,x);
150:   MatMatSolve(G,B,X);
151:   MatDestroy(&G);

153:   /* Out-place Cholesky */
154:   MatGetFactor(Aher,MATSOLVERSCALAPACK,MAT_FACTOR_CHOLESKY,&G);
155:   MatCholeskyFactorSymbolic(G,Aher,0,NULL);
156:   MatCholeskyFactorNumeric(G,Aher,NULL);
157:   if (mats_view) {
158:     MatView(G,PETSC_VIEWER_STDOUT_WORLD);
159:   }
160:   MatSolve(G,b,x);
161:   MatMatSolve(G,B,X);
162:   MatDestroy(&G);

164:   /* Check norm(Aher*x - b) */
165:   VecCreate(PETSC_COMM_WORLD,&c);
166:   VecSetSizes(c,m,PETSC_DECIDE);
167:   VecSetFromOptions(c);
168:   MatMult(Aher,x,c);
169:   VecAXPY(c,-1.0,b);
170:   VecNorm(c,NORM_1,&norm);
171:   if (norm > tol) {
172:     PetscPrintf(PETSC_COMM_WORLD,"Warning: ||Aher*x - b||=%g for Cholesky\n",(double)norm);
173:   }

175:   /* Check norm(Aher*X - B) */
176:   MatMatMult(Aher,X,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C);
177:   MatAXPY(C,-1.0,B,SAME_NONZERO_PATTERN);
178:   MatNorm(C,NORM_1,&norm);
179:   if (norm > tol) {
180:     PetscPrintf(PETSC_COMM_WORLD,"Warning: ||Aher*X - B||=%g for Cholesky\n",(double)norm);
181:   }

183:   /* LU factorization */
184:   /*------------------*/
185:   PetscPrintf(PETSC_COMM_WORLD," Test LU Solver \n");
186:   /* In-place LU */
187:   /* Create matrix factor F, with a copy of A */
188:   MatDuplicate(A,MAT_COPY_VALUES,&F);
189:   /* Create vector d to test MatSolveAdd() */
190:   VecDuplicate(x,&d);
191:   VecCopy(x,d);

193:   /* PF=LU factorization */
194:   MatLUFactor(F,0,0,NULL);

196:   /* Solve LUX = PB */
197:   MatSolveAdd(F,b,d,x);
198:   MatMatSolve(F,B,X);
199:   MatDestroy(&F);

201:   /* Check norm(A*X - B) */
202:   VecCreate(PETSC_COMM_WORLD,&e);
203:   VecSetSizes(e,m,PETSC_DECIDE);
204:   VecSetFromOptions(e);
205:   MatMult(A,x,c);
206:   MatMult(A,d,e);
207:   VecAXPY(c,-1.0,e);
208:   VecAXPY(c,-1.0,b);
209:   VecNorm(c,NORM_1,&norm);
210:   if (norm > tol) {
211:     PetscPrintf(PETSC_COMM_WORLD,"Warning: ||A*x - b||=%g for LU\n",(double)norm);
212:   }
213:   /* Reuse product C; replace Aher with A */
214:   MatProductReplaceMats(A,NULL,NULL,C);
215:   MatMatMult(A,X,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C);
216:   MatAXPY(C,-1.0,B,SAME_NONZERO_PATTERN);
217:   MatNorm(C,NORM_1,&norm);
218:   if (norm > tol) {
219:     PetscPrintf(PETSC_COMM_WORLD,"Warning: ||A*X - B||=%g for LU\n",(double)norm);
220:   }

222:   /* Out-place LU */
223:   MatGetFactor(A,MATSOLVERSCALAPACK,MAT_FACTOR_LU,&F);
224:   MatLUFactorSymbolic(F,A,0,0,NULL);
225:   MatLUFactorNumeric(F,A,NULL);
226:   if (mats_view) {
227:     MatView(F,PETSC_VIEWER_STDOUT_WORLD);
228:   }
229:   MatSolve(F,b,x);
230:   MatMatSolve(F,B,X);
231:   MatDestroy(&F);

233:   /* Free space */
234:   MatDestroy(&A);
235:   MatDestroy(&Aher);
236:   MatDestroy(&B);
237:   MatDestroy(&C);
238:   MatDestroy(&X);
239:   VecDestroy(&b);
240:   VecDestroy(&c);
241:   VecDestroy(&d);
242:   VecDestroy(&e);
243:   VecDestroy(&x);
244:   PetscRandomDestroy(&rand);
245:   PetscFinalize();
246:   return ierr;
247: }

249: /*TEST

251:    build:
252:       requires: scalapack

254:    test:
255:       nsize: 2
256:       output_file: output/ex245.out

258:    test:
259:       suffix: 2
260:       nsize: 6
261:       output_file: output/ex245.out

263: TEST*/