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;
 10:   PetscInt       m=5,n,p,i,j,nrows,ncols;
 11:   PetscScalar    *v,*barray,rval;
 12:   PetscReal      norm,tol=1.e5*PETSC_MACHINE_EPSILON;
 13:   PetscMPIInt    size,rank;
 14:   PetscRandom    rand;
 15:   const PetscInt *rows,*cols;
 16:   IS             isrows,iscols;
 17:   PetscBool      mats_view=PETSC_FALSE;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

248: /*TEST

250:    build:
251:       requires: scalapack

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

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

262: TEST*/