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 { 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]; 17 18 PetscFunctionBeginUser; 19 PetscCall(PetscInitialize(&argc, &args, (char *)0, 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)); 29 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 } 60 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); 72 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)); 81 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)); 87 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 */ 94 95 /* make a symmetric matrix */ 96 if (symm) { 97 Mat AT; 98 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; 107 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")); 115 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]); 133 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)); 165 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)); 170 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 } 184 185 /*TEST 186 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 203 204 TEST*/ 205