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