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 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 PetscCheckFalse(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) { 169 PetscCall(PetscPrintf(PETSC_COMM_SELF,"MatMatSolve: Norm of error %g\n",(double)norm)); 170 } 171 172 /* Free data structures */ 173 PetscCall(MatDestroy(&A)); 174 PetscCall(MatDestroy(&C)); 175 PetscCall(MatDestroy(&F)); 176 PetscCall(MatDestroy(&X)); 177 PetscCall(MatDestroy(&RHS)); 178 PetscCall(PetscRandomDestroy(&rand)); 179 PetscCall(VecDestroy(&x)); 180 PetscCall(VecDestroy(&b)); 181 PetscCall(VecDestroy(&u)); 182 PetscCall(PetscFinalize()); 183 return 0; 184 } 185 186 /*TEST 187 188 testset: 189 output_file: output/ex215.out 190 test: 191 suffix: ns 192 test: 193 suffix: sym 194 args: -symmetric_solve 195 test: 196 suffix: herm 197 args: -hermitian_solve 198 test: 199 suffix: hpd 200 args: -hpd_solve 201 test: 202 suffix: qr 203 args: -ftype qr 204 205 TEST*/ 206