1 2 static char help[] = "Tests cholesky, icc factorization and solve on sequential aij, baij and sbaij matrices. \n"; 3 4 #include <petscmat.h> 5 6 int main(int argc,char **args) 7 { 8 Vec x,y,b; 9 Mat A; /* linear system matrix */ 10 Mat sA,sC; /* symmetric part of the matrices */ 11 PetscInt n,mbs=16,bs=1,nz=3,prob=1,i,j,col[3],block, row,Ii,J,n1,lvl; 12 PetscMPIInt size; 13 PetscReal norm2; 14 PetscScalar neg_one = -1.0,four=4.0,value[3]; 15 IS perm,cperm; 16 PetscRandom rdm; 17 PetscBool reorder = PETSC_FALSE,displ = PETSC_FALSE; 18 MatFactorInfo factinfo; 19 PetscBool equal; 20 PetscBool TestAIJ = PETSC_FALSE,TestBAIJ = PETSC_TRUE; 21 PetscInt TestShift=0; 22 23 PetscFunctionBeginUser; 24 PetscCall(PetscInitialize(&argc,&args,(char*)0,help)); 25 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 26 PetscCheck(size == 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only!"); 27 PetscCall(PetscOptionsGetInt(NULL,NULL,"-bs",&bs,NULL)); 28 PetscCall(PetscOptionsGetInt(NULL,NULL,"-mbs",&mbs,NULL)); 29 PetscCall(PetscOptionsGetBool(NULL,NULL,"-reorder",&reorder,NULL)); 30 PetscCall(PetscOptionsGetBool(NULL,NULL,"-testaij",&TestAIJ,NULL)); 31 PetscCall(PetscOptionsGetInt(NULL,NULL,"-testShift",&TestShift,NULL)); 32 PetscCall(PetscOptionsGetBool(NULL,NULL,"-displ",&displ,NULL)); 33 34 n = mbs*bs; 35 if (TestAIJ) { /* A is in aij format */ 36 PetscCall(MatCreateSeqAIJ(PETSC_COMM_WORLD,n,n,nz,NULL,&A)); 37 TestBAIJ = PETSC_FALSE; 38 } else { /* A is in baij format */ 39 PetscCall(MatCreateSeqBAIJ(PETSC_COMM_WORLD,bs,n,n,nz,NULL,&A)); 40 TestAIJ = PETSC_FALSE; 41 } 42 43 /* Assemble matrix */ 44 if (bs == 1) { 45 PetscCall(PetscOptionsGetInt(NULL,NULL,"-test_problem",&prob,NULL)); 46 if (prob == 1) { /* tridiagonal matrix */ 47 value[0] = -1.0; value[1] = 2.0; value[2] = -1.0; 48 for (i=1; i<n-1; i++) { 49 col[0] = i-1; col[1] = i; col[2] = i+1; 50 PetscCall(MatSetValues(A,1,&i,3,col,value,INSERT_VALUES)); 51 } 52 i = n - 1; col[0]=0; col[1] = n - 2; col[2] = n - 1; 53 54 value[0]= 0.1; value[1]=-1; value[2]=2; 55 PetscCall(MatSetValues(A,1,&i,3,col,value,INSERT_VALUES)); 56 57 i = 0; col[0] = 0; col[1] = 1; col[2]=n-1; 58 59 value[0] = 2.0; value[1] = -1.0; value[2]=0.1; 60 PetscCall(MatSetValues(A,1,&i,3,col,value,INSERT_VALUES)); 61 } else if (prob ==2) { /* matrix for the five point stencil */ 62 n1 = (PetscInt) (PetscSqrtReal((PetscReal)n) + 0.001); 63 PetscCheck(n1*n1 == n,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"sqrt(n) must be a positive integer!"); 64 for (i=0; i<n1; i++) { 65 for (j=0; j<n1; j++) { 66 Ii = j + n1*i; 67 if (i>0) { 68 J = Ii - n1; 69 PetscCall(MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES)); 70 } 71 if (i<n1-1) { 72 J = Ii + n1; 73 PetscCall(MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES)); 74 } 75 if (j>0) { 76 J = Ii - 1; 77 PetscCall(MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES)); 78 } 79 if (j<n1-1) { 80 J = Ii + 1; 81 PetscCall(MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES)); 82 } 83 PetscCall(MatSetValues(A,1,&Ii,1,&Ii,&four,INSERT_VALUES)); 84 } 85 } 86 } 87 } else { /* bs > 1 */ 88 for (block=0; block<n/bs; block++) { 89 /* diagonal blocks */ 90 value[0] = -1.0; value[1] = 4.0; value[2] = -1.0; 91 for (i=1+block*bs; i<bs-1+block*bs; i++) { 92 col[0] = i-1; col[1] = i; col[2] = i+1; 93 PetscCall(MatSetValues(A,1,&i,3,col,value,INSERT_VALUES)); 94 } 95 i = bs - 1+block*bs; col[0] = bs - 2+block*bs; col[1] = bs - 1+block*bs; 96 97 value[0]=-1.0; value[1]=4.0; 98 PetscCall(MatSetValues(A,1,&i,2,col,value,INSERT_VALUES)); 99 100 i = 0+block*bs; col[0] = 0+block*bs; col[1] = 1+block*bs; 101 102 value[0]=4.0; value[1] = -1.0; 103 PetscCall(MatSetValues(A,1,&i,2,col,value,INSERT_VALUES)); 104 } 105 /* off-diagonal blocks */ 106 value[0]=-1.0; 107 for (i=0; i<(n/bs-1)*bs; i++) { 108 col[0]=i+bs; 109 PetscCall(MatSetValues(A,1,&i,1,col,value,INSERT_VALUES)); 110 col[0]=i; row=i+bs; 111 PetscCall(MatSetValues(A,1,&row,1,col,value,INSERT_VALUES)); 112 } 113 } 114 115 if (TestShift) { 116 /* set diagonals in the 0-th block as 0 for testing shift numerical factor */ 117 for (i=0; i<bs; i++) { 118 row = i; col[0] = i; value[0] = 0.0; 119 PetscCall(MatSetValues(A,1,&row,1,col,value,INSERT_VALUES)); 120 } 121 } 122 123 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 124 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 125 126 /* Test MatConvert */ 127 PetscCall(MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE)); 128 PetscCall(MatConvert(A,MATSEQSBAIJ,MAT_INITIAL_MATRIX,&sA)); 129 PetscCall(MatMultEqual(A,sA,20,&equal)); 130 PetscCheck(equal,PETSC_COMM_SELF,PETSC_ERR_USER,"A != sA"); 131 132 /* Test MatGetOwnershipRange() */ 133 PetscCall(MatGetOwnershipRange(A,&Ii,&J)); 134 PetscCall(MatGetOwnershipRange(sA,&i,&j)); 135 PetscCheck(i == Ii && j == J ,PETSC_COMM_SELF,PETSC_ERR_PLIB,"MatGetOwnershipRange() in MatSBAIJ format"); 136 137 /* Vectors */ 138 PetscCall(PetscRandomCreate(PETSC_COMM_SELF,&rdm)); 139 PetscCall(PetscRandomSetFromOptions(rdm)); 140 PetscCall(VecCreateSeq(PETSC_COMM_SELF,n,&x)); 141 PetscCall(VecDuplicate(x,&b)); 142 PetscCall(VecDuplicate(x,&y)); 143 PetscCall(VecSetRandom(x,rdm)); 144 145 /* Test MatReordering() - not work on sbaij matrix */ 146 if (reorder) { 147 PetscCall(MatGetOrdering(A,MATORDERINGRCM,&perm,&cperm)); 148 } else { 149 PetscCall(MatGetOrdering(A,MATORDERINGNATURAL,&perm,&cperm)); 150 } 151 PetscCall(ISDestroy(&cperm)); 152 153 /* initialize factinfo */ 154 PetscCall(MatFactorInfoInitialize(&factinfo)); 155 if (TestShift == 1) { 156 factinfo.shifttype = (PetscReal)MAT_SHIFT_NONZERO; 157 factinfo.shiftamount = 0.1; 158 } else if (TestShift == 2) { 159 factinfo.shifttype = (PetscReal)MAT_SHIFT_POSITIVE_DEFINITE; 160 } 161 162 /* Test MatCholeskyFactor(), MatICCFactor() */ 163 /*------------------------------------------*/ 164 /* Test aij matrix A */ 165 if (TestAIJ) { 166 if (displ) { 167 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"AIJ: \n")); 168 } 169 i = 0; 170 for (lvl=-1; lvl<10; lvl++) { 171 if (lvl==-1) { /* Cholesky factor */ 172 factinfo.fill = 5.0; 173 174 PetscCall(MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&sC)); 175 PetscCall(MatCholeskyFactorSymbolic(sC,A,perm,&factinfo)); 176 } else { /* incomplete Cholesky factor */ 177 factinfo.fill = 5.0; 178 factinfo.levels = lvl; 179 180 PetscCall(MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_ICC,&sC)); 181 PetscCall(MatICCFactorSymbolic(sC,A,perm,&factinfo)); 182 } 183 PetscCall(MatCholeskyFactorNumeric(sC,A,&factinfo)); 184 185 PetscCall(MatMult(A,x,b)); 186 PetscCall(MatSolve(sC,b,y)); 187 PetscCall(MatDestroy(&sC)); 188 189 /* Check the residual */ 190 PetscCall(VecAXPY(y,neg_one,x)); 191 PetscCall(VecNorm(y,NORM_2,&norm2)); 192 193 if (displ) { 194 PetscCall(PetscPrintf(PETSC_COMM_WORLD," lvl: %" PetscInt_FMT ", residual: %g\n", lvl,(double)norm2)); 195 } 196 } 197 } 198 199 /* Test baij matrix A */ 200 if (TestBAIJ) { 201 if (displ) { 202 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"BAIJ: \n")); 203 } 204 i = 0; 205 for (lvl=-1; lvl<10; lvl++) { 206 if (lvl==-1) { /* Cholesky factor */ 207 factinfo.fill = 5.0; 208 209 PetscCall(MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&sC)); 210 PetscCall(MatCholeskyFactorSymbolic(sC,A,perm,&factinfo)); 211 } else { /* incomplete Cholesky factor */ 212 factinfo.fill = 5.0; 213 factinfo.levels = lvl; 214 215 PetscCall(MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_ICC,&sC)); 216 PetscCall(MatICCFactorSymbolic(sC,A,perm,&factinfo)); 217 } 218 PetscCall(MatCholeskyFactorNumeric(sC,A,&factinfo)); 219 220 PetscCall(MatMult(A,x,b)); 221 PetscCall(MatSolve(sC,b,y)); 222 PetscCall(MatDestroy(&sC)); 223 224 /* Check the residual */ 225 PetscCall(VecAXPY(y,neg_one,x)); 226 PetscCall(VecNorm(y,NORM_2,&norm2)); 227 if (displ) { 228 PetscCall(PetscPrintf(PETSC_COMM_WORLD," lvl: %" PetscInt_FMT ", residual: %g\n", lvl,(double)norm2)); 229 } 230 } 231 } 232 233 /* Test sbaij matrix sA */ 234 if (displ) { 235 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"SBAIJ: \n")); 236 } 237 i = 0; 238 for (lvl=-1; lvl<10; lvl++) { 239 if (lvl==-1) { /* Cholesky factor */ 240 factinfo.fill = 5.0; 241 242 PetscCall(MatGetFactor(sA,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&sC)); 243 PetscCall(MatCholeskyFactorSymbolic(sC,sA,perm,&factinfo)); 244 } else { /* incomplete Cholesky factor */ 245 factinfo.fill = 5.0; 246 factinfo.levels = lvl; 247 248 PetscCall(MatGetFactor(sA,MATSOLVERPETSC,MAT_FACTOR_ICC,&sC)); 249 PetscCall(MatICCFactorSymbolic(sC,sA,perm,&factinfo)); 250 } 251 PetscCall(MatCholeskyFactorNumeric(sC,sA,&factinfo)); 252 253 if (lvl==0 && bs==1) { /* Test inplace ICC(0) for sbaij sA - does not work for new datastructure */ 254 /* 255 Mat B; 256 PetscCall(MatDuplicate(sA,MAT_COPY_VALUES,&B)); 257 PetscCall(MatICCFactor(B,perm,&factinfo)); 258 PetscCall(MatEqual(sC,B,&equal)); 259 if (!equal) { 260 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"in-place Cholesky factor != out-place Cholesky factor"); 261 } 262 PetscCall(MatDestroy(&B)); 263 */ 264 } 265 266 PetscCall(MatMult(sA,x,b)); 267 PetscCall(MatSolve(sC,b,y)); 268 269 /* Test MatSolves() */ 270 if (bs == 1) { 271 Vecs xx,bb; 272 PetscCall(VecsCreateSeq(PETSC_COMM_SELF,n,4,&xx)); 273 PetscCall(VecsDuplicate(xx,&bb)); 274 PetscCall(MatSolves(sC,bb,xx)); 275 PetscCall(VecsDestroy(xx)); 276 PetscCall(VecsDestroy(bb)); 277 } 278 PetscCall(MatDestroy(&sC)); 279 280 /* Check the residual */ 281 PetscCall(VecAXPY(y,neg_one,x)); 282 PetscCall(VecNorm(y,NORM_2,&norm2)); 283 if (displ) { 284 PetscCall(PetscPrintf(PETSC_COMM_WORLD," lvl: %" PetscInt_FMT ", residual: %g\n", lvl,(double)norm2)); 285 } 286 } 287 288 PetscCall(ISDestroy(&perm)); 289 PetscCall(MatDestroy(&A)); 290 PetscCall(MatDestroy(&sA)); 291 PetscCall(VecDestroy(&x)); 292 PetscCall(VecDestroy(&y)); 293 PetscCall(VecDestroy(&b)); 294 PetscCall(PetscRandomDestroy(&rdm)); 295 296 PetscCall(PetscFinalize()); 297 return 0; 298 } 299 300 /*TEST 301 302 test: 303 args: -bs {{1 2 3 4 5 6 7 8}} 304 305 test: 306 suffix: 3 307 args: -testaij 308 output_file: output/ex76_1.out 309 310 TEST*/ 311