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