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(PetscStrcpy(file,"${PETSC_DIR}/share/petsc/datafiles/matrices/")); 35 if (hpd) { 36 #if defined(PETSC_USE_COMPLEX) 37 PetscCall(PetscStrcat(file,"hpd-complex-")); 38 #else 39 PetscCall(PetscStrcat(file,"spd-real-")); 40 #endif 41 ftyp = MAT_FACTOR_CHOLESKY; 42 } else { 43 #if defined(PETSC_USE_COMPLEX) 44 PetscCall(PetscStrcat(file,"nh-complex-")); 45 #else 46 PetscCall(PetscStrcat(file,"ns-real-")); 47 #endif 48 } 49 #if defined(PETSC_USE_64BIT_INDICES) 50 PetscCall(PetscStrcat(file,"int64-")); 51 #else 52 PetscCall(PetscStrcat(file,"int32-")); 53 #endif 54 #if defined(PETSC_USE_REAL_SINGLE) 55 PetscCall(PetscStrcat(file,"float32")); 56 #else 57 PetscCall(PetscStrcat(file,"float64")); 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) { 170 PetscCall(PetscPrintf(PETSC_COMM_SELF,"MatMatSolve: Norm of error %g\n",(double)norm)); 171 } 172 173 /* Free data structures */ 174 PetscCall(MatDestroy(&A)); 175 PetscCall(MatDestroy(&C)); 176 PetscCall(MatDestroy(&F)); 177 PetscCall(MatDestroy(&X)); 178 PetscCall(MatDestroy(&RHS)); 179 PetscCall(PetscRandomDestroy(&rand)); 180 PetscCall(VecDestroy(&x)); 181 PetscCall(VecDestroy(&b)); 182 PetscCall(VecDestroy(&u)); 183 PetscCall(PetscFinalize()); 184 return 0; 185 } 186 187 /*TEST 188 189 testset: 190 output_file: output/ex215.out 191 test: 192 suffix: ns 193 test: 194 suffix: sym 195 args: -symmetric_solve 196 test: 197 suffix: herm 198 args: -hermitian_solve 199 test: 200 suffix: hpd 201 args: -hpd_solve 202 test: 203 suffix: qr 204 args: -ftype qr 205 206 TEST*/ 207