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 colums */ 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 if (ftyp == MAT_FACTOR_LU) { 122 ierr = MatLUFactor(F,NULL,NULL,NULL);CHKERRQ(ierr); 123 } else { 124 ierr = MatCholeskyFactor(F,NULL,NULL);CHKERRQ(ierr); 125 } 126 127 for (nsolve = 0; nsolve < 2; nsolve++) { 128 ierr = VecSetRandom(x,rand);CHKERRQ(ierr); 129 ierr = VecCopy(x,u);CHKERRQ(ierr); 130 if (nsolve) { 131 ierr = MatMult(A,x,b);CHKERRQ(ierr); 132 ierr = MatSolve(F,b,x);CHKERRQ(ierr); 133 } else { 134 ierr = MatMultTranspose(A,x,b);CHKERRQ(ierr); 135 ierr = MatSolveTranspose(F,b,x);CHKERRQ(ierr); 136 } 137 /* Check the error */ 138 ierr = VecAXPY(u,-1.0,x);CHKERRQ(ierr); /* u <- (-1.0)x + u */ 139 ierr = VecNorm(u,NORM_2,&norm);CHKERRQ(ierr); 140 if (norm > tol) { 141 PetscReal resi; 142 if (nsolve) { 143 ierr = MatMult(A,x,u);CHKERRQ(ierr); /* u = A*x */ 144 } else { 145 ierr = MatMultTranspose(A,x,u);CHKERRQ(ierr); /* u = A*x */ 146 } 147 ierr = VecAXPY(u,-1.0,b);CHKERRQ(ierr); /* u <- (-1.0)b + u */ 148 ierr = VecNorm(u,NORM_2,&resi);CHKERRQ(ierr); 149 if (nsolve) { 150 ierr = PetscPrintf(PETSC_COMM_SELF,"MatSolve error: Norm of error %g, residual %f\n",norm,resi);CHKERRQ(ierr); 151 } else { 152 ierr = PetscPrintf(PETSC_COMM_SELF,"MatSolveTranspose error: Norm of error %g, residual %f\n",norm,resi);CHKERRQ(ierr); 153 } 154 } 155 } 156 ierr = MatMatMult(A,C,MAT_REUSE_MATRIX,2.0,&RHS);CHKERRQ(ierr); 157 ierr = MatMatSolve(F,RHS,X);CHKERRQ(ierr); 158 159 /* Check the error */ 160 ierr = MatAXPY(X,-1.0,C,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 161 ierr = MatNorm(X,NORM_FROBENIUS,&norm);CHKERRQ(ierr); 162 if (norm > tol) { 163 ierr = PetscPrintf(PETSC_COMM_SELF,"MatMatSolve: Norm of error %g\n",norm);CHKERRQ(ierr); 164 } 165 166 /* Free data structures */ 167 ierr = MatDestroy(&A);CHKERRQ(ierr); 168 ierr = MatDestroy(&C);CHKERRQ(ierr); 169 ierr = MatDestroy(&F);CHKERRQ(ierr); 170 ierr = MatDestroy(&X);CHKERRQ(ierr); 171 ierr = MatDestroy(&RHS);CHKERRQ(ierr); 172 ierr = PetscRandomDestroy(&rand);CHKERRQ(ierr); 173 ierr = VecDestroy(&x);CHKERRQ(ierr); 174 ierr = VecDestroy(&b);CHKERRQ(ierr); 175 ierr = VecDestroy(&u);CHKERRQ(ierr); 176 ierr = PetscFinalize(); 177 return ierr; 178 } 179 180 181 /*TEST 182 183 testset: 184 output_file: output/ex215.out 185 test: 186 suffix: ns 187 test: 188 suffix: sym 189 args: -symmetric_solve 190 test: 191 suffix: herm 192 args: -hermitian_solve 193 test: 194 suffix: hpd 195 args: -hpd_solve 196 197 TEST*/ 198