Actual source code: ex267.c

  1: static char help[] = "Test different MatSolve routines with MATTRANSPOSEVIRTUAL.\n\n";

  3: #include <petscmat.h>

  5: PetscErrorCode TestMatrix(const char *test, Mat A, PetscInt nrhs, PetscBool inplace, PetscBool chol)
  6: {
  7:   Mat       F, RHS, X, C1;
  8:   Vec       b, x, y, f;
  9:   IS        perm, iperm;
 10:   PetscInt  n, i;
 11:   PetscReal norm, tol = 1000 * PETSC_MACHINE_EPSILON;
 12:   PetscBool ht;
 13: #if defined(PETSC_USE_COMPLEX)
 14:   PetscScalar v1 = PetscCMPLX(1.0, -0.1), v2 = PetscCMPLX(-1.0, 0.1);
 15: #else
 16:   PetscScalar v1 = 1.0, v2 = -1.0;
 17: #endif

 19:   PetscFunctionBegin;
 20:   PetscCall(PetscObjectTypeCompare((PetscObject)A, MATHERMITIANTRANSPOSEVIRTUAL, &ht));
 21:   PetscCall(MatCreateVecs(A, &f, &b));
 22:   PetscCall(MatCreateVecs(A, &x, &y));
 23:   PetscCall(VecSet(b, v1));
 24:   PetscCall(VecSet(y, v2));

 26:   PetscCall(MatGetOrdering(A, MATORDERINGND, &perm, &iperm));
 27:   if (!inplace) {
 28:     if (!chol) {
 29:       PetscCall(MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_LU, &F));
 30:       PetscCall(MatLUFactorSymbolic(F, A, perm, iperm, NULL));
 31:       PetscCall(MatLUFactorNumeric(F, A, NULL));
 32:     } else { /* Cholesky */
 33:       PetscCall(MatGetFactor(A, MATSOLVERPETSC, MAT_FACTOR_CHOLESKY, &F));
 34:       PetscCall(MatCholeskyFactorSymbolic(F, A, perm, NULL));
 35:       PetscCall(MatCholeskyFactorNumeric(F, A, NULL));
 36:     }
 37:   } else { /* Test inplace factorization */
 38:     PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &F));
 39:     if (!chol) PetscCall(MatLUFactor(F, perm, iperm, NULL));
 40:     else PetscCall(MatCholeskyFactor(F, perm, NULL));
 41:   }

 43:   /* MatSolve */
 44:   PetscCall(MatSolve(F, b, x));
 45:   PetscCall(MatMult(A, x, f));
 46:   PetscCall(VecAXPY(f, -1.0, b));
 47:   PetscCall(VecNorm(f, NORM_2, &norm));
 48:   if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%12s MatSolve              : Error of norm %g\n", test, (double)norm));

 50:   /* MatSolveTranspose */
 51:   if (!ht) {
 52:     PetscCall(MatSolveTranspose(F, b, x));
 53:     PetscCall(MatMultTranspose(A, x, f));
 54:     PetscCall(VecAXPY(f, -1.0, b));
 55:     PetscCall(VecNorm(f, NORM_2, &norm));
 56:     if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%12s MatSolveTranspose     : Error of norm %g\n", test, (double)norm));
 57:   }

 59:   /* MatSolveAdd */
 60:   PetscCall(MatSolveAdd(F, b, y, x));
 61:   PetscCall(MatMult(A, y, f));
 62:   PetscCall(VecScale(f, -1.0));
 63:   PetscCall(MatMultAdd(A, x, f, f));
 64:   PetscCall(VecAXPY(f, -1.0, b));
 65:   PetscCall(VecNorm(f, NORM_2, &norm));
 66:   if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%12s MatSolveAdd           : Error of norm %g\n", test, (double)norm));

 68:   /* MatSolveTransposeAdd */
 69:   if (!ht) {
 70:     PetscCall(MatSolveTransposeAdd(F, b, y, x));
 71:     PetscCall(MatMultTranspose(A, y, f));
 72:     PetscCall(VecScale(f, -1.0));
 73:     PetscCall(MatMultTransposeAdd(A, x, f, f));
 74:     PetscCall(VecAXPY(f, -1.0, b));
 75:     PetscCall(VecNorm(f, NORM_2, &norm));
 76:     if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%12s MatSolveTransposeAdd  : Error of norm %g\n", test, (double)norm));
 77:   }

 79:   /* MatMatSolve */
 80:   PetscCall(MatGetSize(A, &n, NULL));
 81:   PetscCall(MatCreate(PETSC_COMM_WORLD, &RHS));
 82:   PetscCall(MatSetSizes(RHS, PETSC_DECIDE, PETSC_DECIDE, n, nrhs));
 83:   PetscCall(MatSetType(RHS, MATSEQDENSE));
 84:   PetscCall(MatSetUp(RHS));
 85:   for (i = 0; i < nrhs; i++) PetscCall(MatSetValue(RHS, i, i, 1.0, INSERT_VALUES));
 86:   PetscCall(MatAssemblyBegin(RHS, MAT_FINAL_ASSEMBLY));
 87:   PetscCall(MatAssemblyEnd(RHS, MAT_FINAL_ASSEMBLY));
 88:   PetscCall(MatDuplicate(RHS, MAT_DO_NOT_COPY_VALUES, &X));

 90:   if (!ht) {
 91:     PetscCall(MatMatSolve(F, RHS, X));
 92:     PetscCall(MatMatMult(A, X, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &C1));
 93:     PetscCall(MatAXPY(C1, -1.0, RHS, SAME_NONZERO_PATTERN));
 94:     PetscCall(MatNorm(C1, NORM_FROBENIUS, &norm));
 95:     if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%12s MatMatSolve           : Error of norm %g\n", test, (double)norm));
 96:     PetscCall(MatDestroy(&C1));

 98:     PetscCall(MatMatSolveTranspose(F, RHS, X));
 99:     PetscCall(MatTransposeMatMult(A, X, MAT_INITIAL_MATRIX, PETSC_DETERMINE, &C1));
100:     PetscCall(MatAXPY(C1, -1.0, RHS, SAME_NONZERO_PATTERN));
101:     PetscCall(MatNorm(C1, NORM_FROBENIUS, &norm));
102:     if (norm > tol) PetscCall(PetscPrintf(PETSC_COMM_WORLD, "%12s MatMatSolveTranspose  : Error of norm %g\n", test, (double)norm));
103:     PetscCall(MatDestroy(&C1));
104:   }
105:   PetscCall(VecDestroy(&b));
106:   PetscCall(VecDestroy(&x));
107:   PetscCall(VecDestroy(&f));
108:   PetscCall(VecDestroy(&y));
109:   PetscCall(ISDestroy(&perm));
110:   PetscCall(ISDestroy(&iperm));
111:   PetscCall(MatDestroy(&F));
112:   PetscCall(MatDestroy(&RHS));
113:   PetscCall(MatDestroy(&X));
114:   PetscFunctionReturn(PETSC_SUCCESS);
115: }

117: int main(int argc, char **args)
118: {
119:   PetscMPIInt size;
120:   Mat         A, At, Aht;
121:   PetscInt    i, n = 8, nrhs = 2;
122:   PetscBool   aij, inplace = PETSC_FALSE;
123: #if defined(PETSC_USE_COMPLEX)
124:   PetscScalar a = PetscCMPLX(-1.0, 0.5);
125: #else
126:   PetscScalar a = -1.0;
127: #endif

129:   PetscFunctionBeginUser;
130:   PetscCall(PetscInitialize(&argc, &args, NULL, help));
131:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
132:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor example only");
133:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL));
134:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-nrhs", &nrhs, NULL));
135:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-inplace", &inplace, NULL));
136:   PetscCheck(nrhs <= n, PETSC_COMM_WORLD, PETSC_ERR_ARG_SIZ, "Must have nrhs <= n");

138:   /* Hermitian matrix */
139:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
140:   PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n));
141:   PetscCall(MatSetFromOptions(A));
142:   for (i = 0; i < n; i++) {
143:     if (i > 0) PetscCall(MatSetValue(A, i, i - 1, a, INSERT_VALUES));
144:     if (i < n - 1) PetscCall(MatSetValue(A, i, i + 1, PetscConj(a), INSERT_VALUES));
145:     PetscCall(MatSetValue(A, i, i, 2.0, INSERT_VALUES));
146:   }
147:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
148:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));

150:   PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &aij, MATSEQAIJ, MATSEQBAIJ, ""));
151: #if defined(PETSC_USE_COMPLEX)
152:   PetscCall(MatSetOption(A, MAT_HERMITIAN, PETSC_TRUE));
153: #else
154:   PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_TRUE));
155: #endif

157:   PetscCall(MatCreateTranspose(A, &At));
158:   PetscCall(MatCreateHermitianTranspose(A, &Aht));

160:   PetscCall(TestMatrix("LU T", At, nrhs, inplace, PETSC_FALSE));
161:   PetscCall(TestMatrix("LU HT", Aht, nrhs, inplace, PETSC_FALSE));
162:   if (!aij) {
163:     PetscCall(TestMatrix("Chol T", At, nrhs, inplace, PETSC_TRUE));
164:     PetscCall(TestMatrix("Chol HT", Aht, nrhs, inplace, PETSC_TRUE));
165:   }

167:   /* Make the matrix non-Hermitian */
168:   PetscCall(MatSetValue(A, 0, 1, -5.0, INSERT_VALUES));
169:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
170:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
171: #if defined(PETSC_USE_COMPLEX)
172:   PetscCall(MatSetOption(A, MAT_HERMITIAN, PETSC_FALSE));
173: #else
174:   PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_FALSE));
175: #endif

177:   PetscCall(TestMatrix("LU T nonsym", At, nrhs, inplace, PETSC_FALSE));
178:   PetscCall(TestMatrix("LU HT nonsym", Aht, nrhs, inplace, PETSC_FALSE));

180:   PetscCall(MatDestroy(&A));
181:   PetscCall(MatDestroy(&At));
182:   PetscCall(MatDestroy(&Aht));
183:   PetscCall(PetscFinalize());
184:   return 0;
185: }

187: /*TEST

189:    test:
190:       suffix: 1
191:       args: -inplace {{0 1}} -mat_type {{aij dense}}
192:       output_file: output/empty.out

194: TEST*/