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 CHKERRQ(PetscInitialize(&argc,&args,(char*)0,help)); 19 CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 20 PetscCheckFalse(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 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-symmetric_solve",&symm,NULL)); 26 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-hermitian_solve",&herm,NULL)); 27 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-hpd_solve",&hpd,NULL)); 28 29 /* Determine file from which we read the matrix A */ 30 ftyp = MAT_FACTOR_LU; 31 CHKERRQ(PetscOptionsGetString(NULL,NULL,"-f",file,sizeof(file),&data_provided)); 32 if (!data_provided) { /* get matrices from PETSc distribution */ 33 CHKERRQ(PetscStrcpy(file,"${PETSC_DIR}/share/petsc/datafiles/matrices/")); 34 if (hpd) { 35 #if defined(PETSC_USE_COMPLEX) 36 CHKERRQ(PetscStrcat(file,"hpd-complex-")); 37 #else 38 CHKERRQ(PetscStrcat(file,"spd-real-")); 39 #endif 40 ftyp = MAT_FACTOR_CHOLESKY; 41 } else { 42 #if defined(PETSC_USE_COMPLEX) 43 CHKERRQ(PetscStrcat(file,"nh-complex-")); 44 #else 45 CHKERRQ(PetscStrcat(file,"ns-real-")); 46 #endif 47 } 48 #if defined(PETSC_USE_64BIT_INDICES) 49 CHKERRQ(PetscStrcat(file,"int64-")); 50 #else 51 CHKERRQ(PetscStrcat(file,"int32-")); 52 #endif 53 #if defined(PETSC_USE_REAL_SINGLE) 54 CHKERRQ(PetscStrcat(file,"float32")); 55 #else 56 CHKERRQ(PetscStrcat(file,"float64")); 57 #endif 58 } 59 60 /* Load matrix A */ 61 #if defined(PETSC_USE_REAL___FLOAT128) 62 CHKERRQ(PetscOptionsInsertString(NULL,"-binary_read_double")); 63 #endif 64 CHKERRQ(PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd)); 65 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&A)); 66 CHKERRQ(MatLoad(A,fd)); 67 CHKERRQ(PetscViewerDestroy(&fd)); 68 CHKERRQ(MatConvert(A,MATSEQDENSE,MAT_INPLACE_MATRIX,&A)); 69 CHKERRQ(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 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-nrhs",&nrhs,NULL)); 75 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&C)); 76 CHKERRQ(MatSetSizes(C,m,PETSC_DECIDE,PETSC_DECIDE,nrhs)); 77 CHKERRQ(MatSetType(C,MATDENSE)); 78 CHKERRQ(MatSetFromOptions(C)); 79 CHKERRQ(MatSetUp(C)); 80 81 CHKERRQ(PetscRandomCreate(PETSC_COMM_WORLD,&rand)); 82 CHKERRQ(PetscRandomSetFromOptions(rand)); 83 CHKERRQ(MatSetRandom(C,rand)); 84 CHKERRQ(MatDuplicate(C,MAT_DO_NOT_COPY_VALUES,&X)); 85 CHKERRQ(MatDuplicate(C,MAT_DO_NOT_COPY_VALUES,&RHS)); 86 87 /* Create vectors */ 88 CHKERRQ(VecCreate(PETSC_COMM_WORLD,&x)); 89 CHKERRQ(VecSetSizes(x,n,PETSC_DECIDE)); 90 CHKERRQ(VecSetFromOptions(x)); 91 CHKERRQ(VecDuplicate(x,&b)); 92 CHKERRQ(VecDuplicate(x,&u)); /* save the true solution */ 93 94 /* make a symmetric matrix */ 95 if (symm) { 96 Mat AT; 97 98 CHKERRQ(MatTranspose(A,MAT_INITIAL_MATRIX,&AT)); 99 CHKERRQ(MatAXPY(A,1.0,AT,SAME_NONZERO_PATTERN)); 100 CHKERRQ(MatDestroy(&AT)); 101 ftyp = MAT_FACTOR_CHOLESKY; 102 } 103 /* make an hermitian matrix */ 104 if (herm) { 105 Mat AH; 106 107 CHKERRQ(MatHermitianTranspose(A,MAT_INITIAL_MATRIX,&AH)); 108 CHKERRQ(MatAXPY(A,1.0,AH,SAME_NONZERO_PATTERN)); 109 CHKERRQ(MatDestroy(&AH)); 110 ftyp = MAT_FACTOR_CHOLESKY; 111 } 112 CHKERRQ(PetscObjectSetName((PetscObject)A,"A")); 113 CHKERRQ(MatViewFromOptions(A,NULL,"-amat_view")); 114 115 CHKERRQ(MatDuplicate(A,MAT_COPY_VALUES,&F)); 116 CHKERRQ(MatSetOption(F,MAT_SYMMETRIC,symm)); 117 /* it seems that the SPD concept in PETSc extends naturally to Hermitian Positive definitess */ 118 CHKERRQ(MatSetOption(F,MAT_HERMITIAN,(PetscBool)(hpd || herm))); 119 CHKERRQ(MatSetOption(F,MAT_SPD,hpd)); 120 { 121 PetscInt iftyp = ftyp; 122 CHKERRQ(PetscOptionsGetEList(NULL,NULL,"-ftype",MatFactorTypes,MAT_FACTOR_NUM_TYPES,&iftyp,NULL)); 123 ftyp = (MatFactorType) iftyp; 124 } 125 if (ftyp == MAT_FACTOR_LU) { 126 CHKERRQ(MatLUFactor(F,NULL,NULL,NULL)); 127 } else if (ftyp == MAT_FACTOR_CHOLESKY) { 128 CHKERRQ(MatCholeskyFactor(F,NULL,NULL)); 129 } else if (ftyp == MAT_FACTOR_QR) { 130 CHKERRQ(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 CHKERRQ(VecSetRandom(x,rand)); 135 CHKERRQ(VecCopy(x,u)); 136 if (nsolve) { 137 CHKERRQ(MatMult(A,x,b)); 138 CHKERRQ(MatSolve(F,b,x)); 139 } else { 140 CHKERRQ(MatMultTranspose(A,x,b)); 141 CHKERRQ(MatSolveTranspose(F,b,x)); 142 } 143 /* Check the error */ 144 CHKERRQ(VecAXPY(u,-1.0,x)); /* u <- (-1.0)x + u */ 145 CHKERRQ(VecNorm(u,NORM_2,&norm)); 146 if (norm > tol) { 147 PetscReal resi; 148 if (nsolve) { 149 CHKERRQ(MatMult(A,x,u)); /* u = A*x */ 150 } else { 151 CHKERRQ(MatMultTranspose(A,x,u)); /* u = A*x */ 152 } 153 CHKERRQ(VecAXPY(u,-1.0,b)); /* u <- (-1.0)b + u */ 154 CHKERRQ(VecNorm(u,NORM_2,&resi)); 155 if (nsolve) { 156 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"MatSolve error: Norm of error %g, residual %g\n",(double)norm,(double)resi)); 157 } else { 158 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"MatSolveTranspose error: Norm of error %g, residual %g\n",(double)norm,(double)resi)); 159 } 160 } 161 } 162 CHKERRQ(MatMatMult(A,C,MAT_REUSE_MATRIX,2.0,&RHS)); 163 CHKERRQ(MatMatSolve(F,RHS,X)); 164 165 /* Check the error */ 166 CHKERRQ(MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN)); 167 CHKERRQ(MatNorm(X,NORM_FROBENIUS,&norm)); 168 if (norm > tol) { 169 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"MatMatSolve: Norm of error %g\n",(double)norm)); 170 } 171 172 /* Free data structures */ 173 CHKERRQ(MatDestroy(&A)); 174 CHKERRQ(MatDestroy(&C)); 175 CHKERRQ(MatDestroy(&F)); 176 CHKERRQ(MatDestroy(&X)); 177 CHKERRQ(MatDestroy(&RHS)); 178 CHKERRQ(PetscRandomDestroy(&rand)); 179 CHKERRQ(VecDestroy(&x)); 180 CHKERRQ(VecDestroy(&b)); 181 CHKERRQ(VecDestroy(&u)); 182 CHKERRQ(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