Actual source code: ex245.c

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

  3: #include <petscmat.h>

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

 18:   PetscFunctionBeginUser;
 19:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
 20:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
 21:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

236: /*TEST

238:    build:
239:       requires: scalapack

241:    test:
242:       nsize: 2
243:       requires: !single # garbage prints in single precision from sgemr2d
244:       output_file: output/ex245.out

246:    test:
247:       suffix: 2
248:       nsize: 6
249:       requires: !single # garbage prints in single precision from sgemr2d
250:       output_file: output/ex245.out

252: TEST*/