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