Actual source code: ex215.c

  1: static char help[] = "Tests MatSolve(), MatSolveTranspose() and MatMatSolve() with SEQDENSE\n";

  3: #include <petscmat.h>

  5: int main(int argc, char **args)
  6: {
  7:   Mat           A, RHS, C, F, X;
  8:   Vec           u, x, b;
  9:   PetscMPIInt   size;
 10:   PetscInt      m, n, nsolve, nrhs;
 11:   PetscReal     norm, tol = PETSC_SQRT_MACHINE_EPSILON;
 12:   PetscRandom   rand;
 13:   PetscBool     data_provided, herm, symm, hpd;
 14:   MatFactorType ftyp;
 15:   PetscViewer   fd;
 16:   char          file[PETSC_MAX_PATH_LEN];

 18:   PetscFunctionBeginUser;
 19:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
 20:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
 21:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor test");
 22:   /* Determine which type of solver we want to test for */
 23:   herm = PETSC_FALSE;
 24:   symm = PETSC_FALSE;
 25:   hpd  = PETSC_FALSE;
 26:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-symmetric_solve", &symm, NULL));
 27:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-hermitian_solve", &herm, NULL));
 28:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-hpd_solve", &hpd, NULL));

 30:   /* Determine file from which we read the matrix A */
 31:   ftyp = MAT_FACTOR_LU;
 32:   PetscCall(PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), &data_provided));
 33:   if (!data_provided) { /* get matrices from PETSc distribution */
 34:     PetscCall(PetscStrncpy(file, "${PETSC_DIR}/share/petsc/datafiles/matrices/", sizeof(file)));
 35:     if (hpd) {
 36: #if defined(PETSC_USE_COMPLEX)
 37:       PetscCall(PetscStrlcat(file, "hpd-complex-", sizeof(file)));
 38: #else
 39:       PetscCall(PetscStrlcat(file, "spd-real-", sizeof(file)));
 40: #endif
 41:       ftyp = MAT_FACTOR_CHOLESKY;
 42:     } else {
 43: #if defined(PETSC_USE_COMPLEX)
 44:       PetscCall(PetscStrlcat(file, "nh-complex-", sizeof(file)));
 45: #else
 46:       PetscCall(PetscStrlcat(file, "ns-real-", sizeof(file)));
 47: #endif
 48:     }
 49: #if defined(PETSC_USE_64BIT_INDICES)
 50:     PetscCall(PetscStrlcat(file, "int64-", sizeof(file)));
 51: #else
 52:     PetscCall(PetscStrlcat(file, "int32-", sizeof(file)));
 53: #endif
 54: #if defined(PETSC_USE_REAL_SINGLE)
 55:     PetscCall(PetscStrlcat(file, "float32", sizeof(file)));
 56: #else
 57:     PetscCall(PetscStrlcat(file, "float64", sizeof(file)));
 58: #endif
 59:   }

 61:   /* Load matrix A */
 62: #if defined(PETSC_USE_REAL___FLOAT128)
 63:   PetscCall(PetscOptionsInsertString(NULL, "-binary_read_double"));
 64: #endif
 65:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd));
 66:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
 67:   PetscCall(MatLoad(A, fd));
 68:   PetscCall(PetscViewerDestroy(&fd));
 69:   PetscCall(MatConvert(A, MATSEQDENSE, MAT_INPLACE_MATRIX, &A));
 70:   PetscCall(MatGetSize(A, &m, &n));
 71:   PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%" PetscInt_FMT ", %" PetscInt_FMT ")", m, n);

 73:   /* Create dense matrix C and X; C holds true solution with identical columns */
 74:   nrhs = 2;
 75:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-nrhs", &nrhs, NULL));
 76:   PetscCall(MatCreate(PETSC_COMM_WORLD, &C));
 77:   PetscCall(MatSetSizes(C, m, PETSC_DECIDE, PETSC_DECIDE, nrhs));
 78:   PetscCall(MatSetType(C, MATDENSE));
 79:   PetscCall(MatSetFromOptions(C));
 80:   PetscCall(MatSetUp(C));

 82:   PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand));
 83:   PetscCall(PetscRandomSetFromOptions(rand));
 84:   PetscCall(MatSetRandom(C, rand));
 85:   PetscCall(MatDuplicate(C, MAT_DO_NOT_COPY_VALUES, &X));
 86:   PetscCall(MatDuplicate(C, MAT_DO_NOT_COPY_VALUES, &RHS));

 88:   /* Create vectors */
 89:   PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
 90:   PetscCall(VecSetSizes(x, n, PETSC_DECIDE));
 91:   PetscCall(VecSetFromOptions(x));
 92:   PetscCall(VecDuplicate(x, &b));
 93:   PetscCall(VecDuplicate(x, &u)); /* save the true solution */

 95:   /* make a symmetric matrix */
 96:   if (symm) {
 97:     Mat AT;

 99:     PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &AT));
100:     PetscCall(MatAXPY(A, 1.0, AT, SAME_NONZERO_PATTERN));
101:     PetscCall(MatDestroy(&AT));
102:     ftyp = MAT_FACTOR_CHOLESKY;
103:   }
104:   /* make an hermitian matrix */
105:   if (herm) {
106:     Mat AH;

108:     PetscCall(MatHermitianTranspose(A, MAT_INITIAL_MATRIX, &AH));
109:     PetscCall(MatAXPY(A, 1.0, AH, SAME_NONZERO_PATTERN));
110:     PetscCall(MatDestroy(&AH));
111:     ftyp = MAT_FACTOR_CHOLESKY;
112:   }
113:   PetscCall(PetscObjectSetName((PetscObject)A, "A"));
114:   PetscCall(MatViewFromOptions(A, NULL, "-amat_view"));

116:   PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &F));
117:   PetscCall(MatSetOption(F, MAT_SYMMETRIC, symm));
118:   /* it seems that the SPD concept in PETSc extends naturally to Hermitian Positive definitess */
119:   PetscCall(MatSetOption(F, MAT_HERMITIAN, (PetscBool)(hpd || herm)));
120:   PetscCall(MatSetOption(F, MAT_SPD, hpd));
121:   {
122:     PetscInt iftyp = ftyp;
123:     PetscCall(PetscOptionsGetEList(NULL, NULL, "-ftype", MatFactorTypes, MAT_FACTOR_NUM_TYPES, &iftyp, NULL));
124:     ftyp = (MatFactorType)iftyp;
125:   }
126:   if (ftyp == MAT_FACTOR_LU) {
127:     PetscCall(MatLUFactor(F, NULL, NULL, NULL));
128:   } else if (ftyp == MAT_FACTOR_CHOLESKY) {
129:     PetscCall(MatCholeskyFactor(F, NULL, NULL));
130:   } else if (ftyp == MAT_FACTOR_QR) {
131:     PetscCall(MatQRFactor(F, NULL, NULL));
132:   } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Factorization %s not supported in this example", MatFactorTypes[ftyp]);

134:   for (nsolve = 0; nsolve < 2; nsolve++) {
135:     PetscCall(VecSetRandom(x, rand));
136:     PetscCall(VecCopy(x, u));
137:     if (nsolve) {
138:       PetscCall(MatMult(A, x, b));
139:       PetscCall(MatSolve(F, b, x));
140:     } else {
141:       PetscCall(MatMultTranspose(A, x, b));
142:       PetscCall(MatSolveTranspose(F, b, x));
143:     }
144:     /* Check the error */
145:     PetscCall(VecAXPY(u, -1.0, x)); /* u <- (-1.0)x + u */
146:     PetscCall(VecNorm(u, NORM_2, &norm));
147:     if (norm > tol) {
148:       PetscReal resi;
149:       if (nsolve) {
150:         PetscCall(MatMult(A, x, u)); /* u = A*x */
151:       } else {
152:         PetscCall(MatMultTranspose(A, x, u)); /* u = A*x */
153:       }
154:       PetscCall(VecAXPY(u, -1.0, b)); /* u <- (-1.0)b + u */
155:       PetscCall(VecNorm(u, NORM_2, &resi));
156:       if (nsolve) {
157:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "MatSolve error: Norm of error %g, residual %g\n", (double)norm, (double)resi));
158:       } else {
159:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "MatSolveTranspose error: Norm of error %g, residual %g\n", (double)norm, (double)resi));
160:       }
161:     }
162:   }
163:   PetscCall(MatMatMult(A, C, MAT_REUSE_MATRIX, 2.0, &RHS));
164:   PetscCall(MatMatSolve(F, RHS, X));

166:   /* Check the error */
167:   PetscCall(MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN));
168:   PetscCall(MatNorm(X, NORM_FROBENIUS, &norm));
169:   if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_SELF, "MatMatSolve: Norm of error %g\n", (double)norm));

171:   /* Free data structures */
172:   PetscCall(MatDestroy(&A));
173:   PetscCall(MatDestroy(&C));
174:   PetscCall(MatDestroy(&F));
175:   PetscCall(MatDestroy(&X));
176:   PetscCall(MatDestroy(&RHS));
177:   PetscCall(PetscRandomDestroy(&rand));
178:   PetscCall(VecDestroy(&x));
179:   PetscCall(VecDestroy(&b));
180:   PetscCall(VecDestroy(&u));
181:   PetscCall(PetscFinalize());
182:   return 0;
183: }

185: /*TEST

187:   testset:
188:     output_file: output/ex215.out
189:     test:
190:       suffix: ns
191:     test:
192:       suffix: sym
193:       args: -symmetric_solve
194:     test:
195:       suffix: herm
196:       args: -hermitian_solve
197:     test:
198:       suffix: hpd
199:       args: -hpd_solve
200:     test:
201:       suffix: qr
202:       args: -ftype qr

204: TEST*/