1 static char help[] = "Tests MatSolve(), MatSolveTranspose() and MatMatSolve() with SEQDENSE\n"; 2 3 #include <petscmat.h> 4 5 int main(int argc, char **args) { 6 Mat A, RHS, C, F, X; 7 Vec u, x, b; 8 PetscMPIInt size; 9 PetscInt m, n, nsolve, nrhs; 10 PetscReal norm, tol = PETSC_SQRT_MACHINE_EPSILON; 11 PetscRandom rand; 12 PetscBool data_provided, herm, symm, hpd; 13 MatFactorType ftyp; 14 PetscViewer fd; 15 char file[PETSC_MAX_PATH_LEN]; 16 17 PetscFunctionBeginUser; 18 PetscCall(PetscInitialize(&argc, &args, (char *)0, help)); 19 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 20 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This is a uniprocessor test"); 21 /* Determine which type of solver we want to test for */ 22 herm = PETSC_FALSE; 23 symm = PETSC_FALSE; 24 hpd = PETSC_FALSE; 25 PetscCall(PetscOptionsGetBool(NULL, NULL, "-symmetric_solve", &symm, NULL)); 26 PetscCall(PetscOptionsGetBool(NULL, NULL, "-hermitian_solve", &herm, NULL)); 27 PetscCall(PetscOptionsGetBool(NULL, NULL, "-hpd_solve", &hpd, NULL)); 28 29 /* Determine file from which we read the matrix A */ 30 ftyp = MAT_FACTOR_LU; 31 PetscCall(PetscOptionsGetString(NULL, NULL, "-f", file, sizeof(file), &data_provided)); 32 if (!data_provided) { /* get matrices from PETSc distribution */ 33 PetscCall(PetscStrcpy(file, "${PETSC_DIR}/share/petsc/datafiles/matrices/")); 34 if (hpd) { 35 #if defined(PETSC_USE_COMPLEX) 36 PetscCall(PetscStrcat(file, "hpd-complex-")); 37 #else 38 PetscCall(PetscStrcat(file, "spd-real-")); 39 #endif 40 ftyp = MAT_FACTOR_CHOLESKY; 41 } else { 42 #if defined(PETSC_USE_COMPLEX) 43 PetscCall(PetscStrcat(file, "nh-complex-")); 44 #else 45 PetscCall(PetscStrcat(file, "ns-real-")); 46 #endif 47 } 48 #if defined(PETSC_USE_64BIT_INDICES) 49 PetscCall(PetscStrcat(file, "int64-")); 50 #else 51 PetscCall(PetscStrcat(file, "int32-")); 52 #endif 53 #if defined(PETSC_USE_REAL_SINGLE) 54 PetscCall(PetscStrcat(file, "float32")); 55 #else 56 PetscCall(PetscStrcat(file, "float64")); 57 #endif 58 } 59 60 /* Load matrix A */ 61 #if defined(PETSC_USE_REAL___FLOAT128) 62 PetscCall(PetscOptionsInsertString(NULL, "-binary_read_double")); 63 #endif 64 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, file, FILE_MODE_READ, &fd)); 65 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 66 PetscCall(MatLoad(A, fd)); 67 PetscCall(PetscViewerDestroy(&fd)); 68 PetscCall(MatConvert(A, MATSEQDENSE, MAT_INPLACE_MATRIX, &A)); 69 PetscCall(MatGetSize(A, &m, &n)); 70 PetscCheck(m == n, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "This example is not intended for rectangular matrices (%" PetscInt_FMT ", %" PetscInt_FMT ")", m, n); 71 72 /* Create dense matrix C and X; C holds true solution with identical columns */ 73 nrhs = 2; 74 PetscCall(PetscOptionsGetInt(NULL, NULL, "-nrhs", &nrhs, NULL)); 75 PetscCall(MatCreate(PETSC_COMM_WORLD, &C)); 76 PetscCall(MatSetSizes(C, m, PETSC_DECIDE, PETSC_DECIDE, nrhs)); 77 PetscCall(MatSetType(C, MATDENSE)); 78 PetscCall(MatSetFromOptions(C)); 79 PetscCall(MatSetUp(C)); 80 81 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rand)); 82 PetscCall(PetscRandomSetFromOptions(rand)); 83 PetscCall(MatSetRandom(C, rand)); 84 PetscCall(MatDuplicate(C, MAT_DO_NOT_COPY_VALUES, &X)); 85 PetscCall(MatDuplicate(C, MAT_DO_NOT_COPY_VALUES, &RHS)); 86 87 /* Create vectors */ 88 PetscCall(VecCreate(PETSC_COMM_WORLD, &x)); 89 PetscCall(VecSetSizes(x, n, PETSC_DECIDE)); 90 PetscCall(VecSetFromOptions(x)); 91 PetscCall(VecDuplicate(x, &b)); 92 PetscCall(VecDuplicate(x, &u)); /* save the true solution */ 93 94 /* make a symmetric matrix */ 95 if (symm) { 96 Mat AT; 97 98 PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &AT)); 99 PetscCall(MatAXPY(A, 1.0, AT, SAME_NONZERO_PATTERN)); 100 PetscCall(MatDestroy(&AT)); 101 ftyp = MAT_FACTOR_CHOLESKY; 102 } 103 /* make an hermitian matrix */ 104 if (herm) { 105 Mat AH; 106 107 PetscCall(MatHermitianTranspose(A, MAT_INITIAL_MATRIX, &AH)); 108 PetscCall(MatAXPY(A, 1.0, AH, SAME_NONZERO_PATTERN)); 109 PetscCall(MatDestroy(&AH)); 110 ftyp = MAT_FACTOR_CHOLESKY; 111 } 112 PetscCall(PetscObjectSetName((PetscObject)A, "A")); 113 PetscCall(MatViewFromOptions(A, NULL, "-amat_view")); 114 115 PetscCall(MatDuplicate(A, MAT_COPY_VALUES, &F)); 116 PetscCall(MatSetOption(F, MAT_SYMMETRIC, symm)); 117 /* it seems that the SPD concept in PETSc extends naturally to Hermitian Positive definitess */ 118 PetscCall(MatSetOption(F, MAT_HERMITIAN, (PetscBool)(hpd || herm))); 119 PetscCall(MatSetOption(F, MAT_SPD, hpd)); 120 { 121 PetscInt iftyp = ftyp; 122 PetscCall(PetscOptionsGetEList(NULL, NULL, "-ftype", MatFactorTypes, MAT_FACTOR_NUM_TYPES, &iftyp, NULL)); 123 ftyp = (MatFactorType)iftyp; 124 } 125 if (ftyp == MAT_FACTOR_LU) { 126 PetscCall(MatLUFactor(F, NULL, NULL, NULL)); 127 } else if (ftyp == MAT_FACTOR_CHOLESKY) { 128 PetscCall(MatCholeskyFactor(F, NULL, NULL)); 129 } else if (ftyp == MAT_FACTOR_QR) { 130 PetscCall(MatQRFactor(F, NULL, NULL)); 131 } else SETERRQ(PETSC_COMM_WORLD, PETSC_ERR_SUP, "Factorization %s not supported in this example", MatFactorTypes[ftyp]); 132 133 for (nsolve = 0; nsolve < 2; nsolve++) { 134 PetscCall(VecSetRandom(x, rand)); 135 PetscCall(VecCopy(x, u)); 136 if (nsolve) { 137 PetscCall(MatMult(A, x, b)); 138 PetscCall(MatSolve(F, b, x)); 139 } else { 140 PetscCall(MatMultTranspose(A, x, b)); 141 PetscCall(MatSolveTranspose(F, b, x)); 142 } 143 /* Check the error */ 144 PetscCall(VecAXPY(u, -1.0, x)); /* u <- (-1.0)x + u */ 145 PetscCall(VecNorm(u, NORM_2, &norm)); 146 if (norm > tol) { 147 PetscReal resi; 148 if (nsolve) { 149 PetscCall(MatMult(A, x, u)); /* u = A*x */ 150 } else { 151 PetscCall(MatMultTranspose(A, x, u)); /* u = A*x */ 152 } 153 PetscCall(VecAXPY(u, -1.0, b)); /* u <- (-1.0)b + u */ 154 PetscCall(VecNorm(u, NORM_2, &resi)); 155 if (nsolve) { 156 PetscCall(PetscPrintf(PETSC_COMM_SELF, "MatSolve error: Norm of error %g, residual %g\n", (double)norm, (double)resi)); 157 } else { 158 PetscCall(PetscPrintf(PETSC_COMM_SELF, "MatSolveTranspose error: Norm of error %g, residual %g\n", (double)norm, (double)resi)); 159 } 160 } 161 } 162 PetscCall(MatMatMult(A, C, MAT_REUSE_MATRIX, 2.0, &RHS)); 163 PetscCall(MatMatSolve(F, RHS, X)); 164 165 /* Check the error */ 166 PetscCall(MatAXPY(X, -1.0, C, SAME_NONZERO_PATTERN)); 167 PetscCall(MatNorm(X, NORM_FROBENIUS, &norm)); 168 if (norm > tol) { PetscCall(PetscPrintf(PETSC_COMM_SELF, "MatMatSolve: Norm of error %g\n", (double)norm)); } 169 170 /* Free data structures */ 171 PetscCall(MatDestroy(&A)); 172 PetscCall(MatDestroy(&C)); 173 PetscCall(MatDestroy(&F)); 174 PetscCall(MatDestroy(&X)); 175 PetscCall(MatDestroy(&RHS)); 176 PetscCall(PetscRandomDestroy(&rand)); 177 PetscCall(VecDestroy(&x)); 178 PetscCall(VecDestroy(&b)); 179 PetscCall(VecDestroy(&u)); 180 PetscCall(PetscFinalize()); 181 return 0; 182 } 183 184 /*TEST 185 186 testset: 187 output_file: output/ex215.out 188 test: 189 suffix: ns 190 test: 191 suffix: sym 192 args: -symmetric_solve 193 test: 194 suffix: herm 195 args: -hermitian_solve 196 test: 197 suffix: hpd 198 args: -hpd_solve 199 test: 200 suffix: qr 201 args: -ftype qr 202 203 TEST*/ 204