1 static char help[] = "Test different MatSolve routines with MATTRANSPOSEVIRTUAL.\n\n"; 2 3 #include <petscmat.h> 4 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 18 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)); 25 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 } 42 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)); 49 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 } 58 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)); 67 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 } 78 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)); 89 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)); 97 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 } 116 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 128 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"); 137 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)); 149 150 PetscCall(PetscObjectTypeCompareAny((PetscObject)A, &aij, MATSEQAIJ, MATSEQBAIJ, "")); 151 PetscCall(MatSetOption(A, PetscDefined(USE_COMPLEX) ? MAT_HERMITIAN : MAT_SYMMETRIC, PETSC_TRUE)); 152 153 PetscCall(MatCreateTranspose(A, &At)); 154 PetscCall(MatCreateHermitianTranspose(A, &Aht)); 155 156 PetscCall(TestMatrix("LU T", At, nrhs, inplace, PETSC_FALSE)); 157 PetscCall(TestMatrix("LU HT", Aht, nrhs, inplace, PETSC_FALSE)); 158 if (!aij) { 159 PetscCall(TestMatrix("Chol T", At, nrhs, inplace, PETSC_TRUE)); 160 PetscCall(TestMatrix("Chol HT", Aht, nrhs, inplace, PETSC_TRUE)); 161 } 162 163 /* Make the matrix non-Hermitian */ 164 PetscCall(MatSetValue(A, 0, 1, -5.0, INSERT_VALUES)); 165 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 166 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 167 #if defined(PETSC_USE_COMPLEX) 168 PetscCall(MatSetOption(A, MAT_HERMITIAN, PETSC_FALSE)); 169 #else 170 PetscCall(MatSetOption(A, MAT_SYMMETRIC, PETSC_FALSE)); 171 #endif 172 173 PetscCall(TestMatrix("LU T nonsym", At, nrhs, inplace, PETSC_FALSE)); 174 PetscCall(TestMatrix("LU HT nonsym", Aht, nrhs, inplace, PETSC_FALSE)); 175 176 PetscCall(MatDestroy(&A)); 177 PetscCall(MatDestroy(&At)); 178 PetscCall(MatDestroy(&Aht)); 179 PetscCall(PetscFinalize()); 180 return 0; 181 } 182 183 /*TEST 184 185 test: 186 suffix: 1 187 args: -inplace {{0 1}} -mat_type {{aij dense}} 188 output_file: output/empty.out 189 190 TEST*/ 191