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];

 19:   PetscInitialize(&argc, &args, (char *)0, help);
 20:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
 22:   /* Determine which type of solver we want to test for */
 23:   herm = PETSC_FALSE;
 24:   symm = PETSC_FALSE;
 25:   hpd  = PETSC_FALSE;
 26:   PetscOptionsGetBool(NULL, NULL, "-symmetric_solve", &symm, NULL);
 27:   PetscOptionsGetBool(NULL, NULL, "-hermitian_solve", &herm, NULL);
 28:   PetscOptionsGetBool(NULL, NULL, "-hpd_solve", &hpd, NULL);

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

 61:   /* Load matrix A */
 62: #if defined(PETSC_USE_REAL___FLOAT128)
 63:   PetscOptionsInsertString(NULL, "-binary_read_double");
 64: #endif
 65:   PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd);
 66:   MatCreate(PETSC_COMM_WORLD, &A);
 67:   MatLoad(A, fd);
 68:   PetscViewerDestroy(&fd);
 69:   MatConvert(A, MATSEQDENSE, MAT_INPLACE_MATRIX, &A);
 70:   MatGetSize(A, &m, &n);

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

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

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

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

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

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

116:   MatDuplicate(A, MAT_COPY_VALUES, &F);
117:   MatSetOption(F, MAT_SYMMETRIC, symm);
118:   /* it seems that the SPD concept in PETSc extends naturally to Hermitian Positive definitess */
119:   MatSetOption(F, MAT_HERMITIAN, (PetscBool)(hpd || herm));
120:   MatSetOption(F, MAT_SPD, hpd);
121:   {
122:     PetscInt iftyp = ftyp;
123:     PetscOptionsGetEList(NULL, NULL, "-ftype", MatFactorTypes, MAT_FACTOR_NUM_TYPES, &iftyp, NULL);
124:     ftyp = (MatFactorType)iftyp;
125:   }
126:   if (ftyp == MAT_FACTOR_LU) {
127:     MatLUFactor(F, NULL, NULL, NULL);
128:   } else if (ftyp == MAT_FACTOR_CHOLESKY) {
129:     MatCholeskyFactor(F, NULL, NULL);
130:   } else if (ftyp == MAT_FACTOR_QR) {
131:     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:     VecSetRandom(x, rand);
136:     VecCopy(x, u);
137:     if (nsolve) {
138:       MatMult(A, x, b);
139:       MatSolve(F, b, x);
140:     } else {
141:       MatMultTranspose(A, x, b);
142:       MatSolveTranspose(F, b, x);
143:     }
144:     /* Check the error */
145:     VecAXPY(u, -1.0, x); /* u <- (-1.0)x + u */
146:     VecNorm(u, NORM_2, &norm);
147:     if (norm > tol) {
148:       PetscReal resi;
149:       if (nsolve) {
150:         MatMult(A, x, u); /* u = A*x */
151:       } else {
152:         MatMultTranspose(A, x, u); /* u = A*x */
153:       }
154:       VecAXPY(u, -1.0, b); /* u <- (-1.0)b + u */
155:       VecNorm(u, NORM_2, &resi);
156:       if (nsolve) {
157:         PetscPrintf(PETSC_COMM_SELF, "MatSolve error: Norm of error %g, residual %g\n", (double)norm, (double)resi);
158:       } else {
159:         PetscPrintf(PETSC_COMM_SELF, "MatSolveTranspose error: Norm of error %g, residual %g\n", (double)norm, (double)resi);
160:       }
161:     }
162:   }
163:   MatMatMult(A, C, MAT_REUSE_MATRIX, 2.0, &RHS);
164:   MatMatSolve(F, RHS, X);

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

171:   /* Free data structures */
172:   MatDestroy(&A);
173:   MatDestroy(&C);
174:   MatDestroy(&F);
175:   MatDestroy(&X);
176:   MatDestroy(&RHS);
177:   PetscRandomDestroy(&rand);
178:   VecDestroy(&x);
179:   VecDestroy(&b);
180:   VecDestroy(&u);
181:   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*/