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:   PetscFunctionBeginUser;
 20:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
 21:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
 22:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));

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

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

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

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

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

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

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

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

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

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

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

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

171:   /* Check norm(Aher*X - B) */
172:   PetscCall(MatMatMult(Aher, X, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &C));
173:   PetscCall(MatAXPY(C, -1.0, B, SAME_NONZERO_PATTERN));
174:   PetscCall(MatNorm(C, NORM_1, &norm));
175:   if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Warning: ||Aher*X - B||=%g for Cholesky\n", (double)norm));

177:   /* LU factorization */
178:   /*------------------*/
179:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, " Test LU Solver \n"));
180:   /* In-place LU */
181:   /* Create matrix factor F, with a copy of A */
182:   PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &F));
183:   /* Create vector d to test MatSolveAdd() */
184:   PetscCall(VecDuplicate(x, &d));
185:   PetscCall(VecCopy(x, d));

187:   /* PF=LU factorization */
188:   PetscCall(MatLUFactor(F, 0, 0, NULL));

190:   /* Solve LUX = PB */
191:   PetscCall(MatSolveAdd(F, b, d, x));
192:   PetscCall(MatMatSolve(F, B, X));
193:   PetscCall(MatDestroy(&F));

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

212:   /* Out-place LU */
213:   PetscCall(MatGetFactor(A, MATSOLVERSCALAPACK, MAT_FACTOR_LU, &F));
214:   PetscCall(MatLUFactorSymbolic(F, A, 0, 0, NULL));
215:   PetscCall(MatLUFactorNumeric(F, A, NULL));
216:   if (mats_view) PetscCall(MatView(F, PETSC_VIEWER_STDOUT_WORLD));
217:   PetscCall(MatSolve(F, b, x));
218:   PetscCall(MatMatSolve(F, B, X));
219:   PetscCall(MatDestroy(&F));

221:   /* Free space */
222:   PetscCall(MatDestroy(&A));
223:   PetscCall(MatDestroy(&Aher));
224:   PetscCall(MatDestroy(&B));
225:   PetscCall(MatDestroy(&C));
226:   PetscCall(MatDestroy(&X));
227:   PetscCall(VecDestroy(&b));
228:   PetscCall(VecDestroy(&c));
229:   PetscCall(VecDestroy(&d));
230:   PetscCall(VecDestroy(&e));
231:   PetscCall(VecDestroy(&x));
232:   PetscCall(PetscRandomDestroy(&rand));
233:   PetscCall(PetscFinalize());
234:   return 0;
235: }

237: /*TEST

239:    build:
240:       requires: scalapack

242:    test:
243:       nsize: 2
244:       output_file: output/ex245.out

246:    test:
247:       suffix: 2
248:       nsize: 6
249:       output_file: output/ex245.out

251: TEST*/