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 PetscErrorCode ierr; 13 PetscMPIInt size; 14 PetscReal norm2; 15 PetscScalar neg_one = -1.0,four=4.0,value[3]; 16 IS perm,cperm; 17 PetscRandom rdm; 18 PetscBool reorder = PETSC_FALSE,displ = PETSC_FALSE; 19 MatFactorInfo factinfo; 20 PetscBool equal; 21 PetscBool TestAIJ = PETSC_FALSE,TestBAIJ = PETSC_TRUE; 22 PetscInt TestShift=0; 23 24 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; 25 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr); 26 if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only!"); 27 ierr = PetscOptionsGetInt(NULL,NULL,"-bs",&bs,NULL);CHKERRQ(ierr); 28 ierr = PetscOptionsGetInt(NULL,NULL,"-mbs",&mbs,NULL);CHKERRQ(ierr); 29 ierr = PetscOptionsGetBool(NULL,NULL,"-reorder",&reorder,NULL);CHKERRQ(ierr); 30 ierr = PetscOptionsGetBool(NULL,NULL,"-testaij",&TestAIJ,NULL);CHKERRQ(ierr); 31 ierr = PetscOptionsGetInt(NULL,NULL,"-testShift",&TestShift,NULL);CHKERRQ(ierr); 32 ierr = PetscOptionsGetBool(NULL,NULL,"-displ",&displ,NULL);CHKERRQ(ierr); 33 34 n = mbs*bs; 35 if (TestAIJ) { /* A is in aij format */ 36 ierr = MatCreateSeqAIJ(PETSC_COMM_WORLD,n,n,nz,NULL,&A);CHKERRQ(ierr); 37 TestBAIJ = PETSC_FALSE; 38 } else { /* A is in baij format */ 39 ierr =MatCreateSeqBAIJ(PETSC_COMM_WORLD,bs,n,n,nz,NULL,&A);CHKERRQ(ierr); 40 TestAIJ = PETSC_FALSE; 41 } 42 43 /* Assemble matrix */ 44 if (bs == 1) { 45 ierr = PetscOptionsGetInt(NULL,NULL,"-test_problem",&prob,NULL);CHKERRQ(ierr); 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 ierr = MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);CHKERRQ(ierr); 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 ierr = MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);CHKERRQ(ierr); 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 ierr = MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);CHKERRQ(ierr); 61 } else if (prob ==2) { /* matrix for the five point stencil */ 62 n1 = (PetscInt) (PetscSqrtReal((PetscReal)n) + 0.001); 63 if (n1*n1 - n) SETERRQ(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 ierr = MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);CHKERRQ(ierr); 70 } 71 if (i<n1-1) { 72 J = Ii + n1; 73 ierr = MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);CHKERRQ(ierr); 74 } 75 if (j>0) { 76 J = Ii - 1; 77 ierr = MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);CHKERRQ(ierr); 78 } 79 if (j<n1-1) { 80 J = Ii + 1; 81 ierr = MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES);CHKERRQ(ierr); 82 } 83 ierr = MatSetValues(A,1,&Ii,1,&Ii,&four,INSERT_VALUES);CHKERRQ(ierr); 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 ierr = MatSetValues(A,1,&i,3,col,value,INSERT_VALUES);CHKERRQ(ierr); 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 ierr = MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);CHKERRQ(ierr); 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 ierr = MatSetValues(A,1,&i,2,col,value,INSERT_VALUES);CHKERRQ(ierr); 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 ierr = MatSetValues(A,1,&i,1,col,value,INSERT_VALUES);CHKERRQ(ierr); 110 col[0]=i; row=i+bs; 111 ierr = MatSetValues(A,1,&row,1,col,value,INSERT_VALUES);CHKERRQ(ierr); 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 ierr = MatSetValues(A,1,&row,1,col,value,INSERT_VALUES);CHKERRQ(ierr); 120 } 121 } 122 123 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 124 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 125 126 /* Test MatConvert */ 127 ierr = MatSetOption(A,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr); 128 ierr = MatConvert(A,MATSEQSBAIJ,MAT_INITIAL_MATRIX,&sA);CHKERRQ(ierr); 129 ierr = MatMultEqual(A,sA,20,&equal);CHKERRQ(ierr); 130 if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"A != sA"); 131 132 /* Test MatGetOwnershipRange() */ 133 ierr = MatGetOwnershipRange(A,&Ii,&J);CHKERRQ(ierr); 134 ierr = MatGetOwnershipRange(sA,&i,&j);CHKERRQ(ierr); 135 if (i-Ii || j-J) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"MatGetOwnershipRange() in MatSBAIJ format\n"); 136 137 /* Vectors */ 138 ierr = PetscRandomCreate(PETSC_COMM_SELF,&rdm);CHKERRQ(ierr); 139 ierr = PetscRandomSetFromOptions(rdm);CHKERRQ(ierr); 140 ierr = VecCreateSeq(PETSC_COMM_SELF,n,&x);CHKERRQ(ierr); 141 ierr = VecDuplicate(x,&b);CHKERRQ(ierr); 142 ierr = VecDuplicate(x,&y);CHKERRQ(ierr); 143 ierr = VecSetRandom(x,rdm);CHKERRQ(ierr); 144 145 /* Test MatReordering() - not work on sbaij matrix */ 146 if (reorder) { 147 ierr = MatGetOrdering(A,MATORDERINGRCM,&perm,&cperm);CHKERRQ(ierr); 148 } else { 149 ierr = MatGetOrdering(A,MATORDERINGNATURAL,&perm,&cperm);CHKERRQ(ierr); 150 } 151 ierr = ISDestroy(&cperm);CHKERRQ(ierr); 152 153 /* initialize factinfo */ 154 ierr = MatFactorInfoInitialize(&factinfo);CHKERRQ(ierr); 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 ierr = PetscPrintf(PETSC_COMM_WORLD,"AIJ: \n");CHKERRQ(ierr); 168 } 169 i = 0; 170 for (lvl=-1; lvl<10; lvl++) { 171 if (lvl==-1) { /* Cholesky factor */ 172 factinfo.fill = 5.0; 173 174 ierr = MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&sC);CHKERRQ(ierr); 175 ierr = MatCholeskyFactorSymbolic(sC,A,perm,&factinfo);CHKERRQ(ierr); 176 } else { /* incomplete Cholesky factor */ 177 factinfo.fill = 5.0; 178 factinfo.levels = lvl; 179 180 ierr = MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_ICC,&sC);CHKERRQ(ierr); 181 ierr = MatICCFactorSymbolic(sC,A,perm,&factinfo);CHKERRQ(ierr); 182 } 183 ierr = MatCholeskyFactorNumeric(sC,A,&factinfo);CHKERRQ(ierr); 184 185 ierr = MatMult(A,x,b);CHKERRQ(ierr); 186 ierr = MatSolve(sC,b,y);CHKERRQ(ierr); 187 ierr = MatDestroy(&sC);CHKERRQ(ierr); 188 189 /* Check the residual */ 190 ierr = VecAXPY(y,neg_one,x);CHKERRQ(ierr); 191 ierr = VecNorm(y,NORM_2,&norm2);CHKERRQ(ierr); 192 193 if (displ) { 194 ierr = PetscPrintf(PETSC_COMM_WORLD," lvl: %D, residual: %g\n", lvl,(double)norm2);CHKERRQ(ierr); 195 } 196 } 197 } 198 199 /* Test baij matrix A */ 200 if (TestBAIJ) { 201 if (displ) { 202 ierr = PetscPrintf(PETSC_COMM_WORLD,"BAIJ: \n");CHKERRQ(ierr); 203 } 204 i = 0; 205 for (lvl=-1; lvl<10; lvl++) { 206 if (lvl==-1) { /* Cholesky factor */ 207 factinfo.fill = 5.0; 208 209 ierr = MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&sC);CHKERRQ(ierr); 210 ierr = MatCholeskyFactorSymbolic(sC,A,perm,&factinfo);CHKERRQ(ierr); 211 } else { /* incomplete Cholesky factor */ 212 factinfo.fill = 5.0; 213 factinfo.levels = lvl; 214 215 ierr = MatGetFactor(A,MATSOLVERPETSC,MAT_FACTOR_ICC,&sC);CHKERRQ(ierr); 216 ierr = MatICCFactorSymbolic(sC,A,perm,&factinfo);CHKERRQ(ierr); 217 } 218 ierr = MatCholeskyFactorNumeric(sC,A,&factinfo);CHKERRQ(ierr); 219 220 ierr = MatMult(A,x,b);CHKERRQ(ierr); 221 ierr = MatSolve(sC,b,y);CHKERRQ(ierr); 222 ierr = MatDestroy(&sC);CHKERRQ(ierr); 223 224 /* Check the residual */ 225 ierr = VecAXPY(y,neg_one,x);CHKERRQ(ierr); 226 ierr = VecNorm(y,NORM_2,&norm2);CHKERRQ(ierr); 227 if (displ) { 228 ierr = PetscPrintf(PETSC_COMM_WORLD," lvl: %D, residual: %g\n", lvl,(double)norm2);CHKERRQ(ierr); 229 } 230 } 231 } 232 233 /* Test sbaij matrix sA */ 234 if (displ) { 235 ierr = PetscPrintf(PETSC_COMM_WORLD,"SBAIJ: \n");CHKERRQ(ierr); 236 } 237 i = 0; 238 for (lvl=-1; lvl<10; lvl++) { 239 if (lvl==-1) { /* Cholesky factor */ 240 factinfo.fill = 5.0; 241 242 ierr = MatGetFactor(sA,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&sC);CHKERRQ(ierr); 243 ierr = MatCholeskyFactorSymbolic(sC,sA,perm,&factinfo);CHKERRQ(ierr); 244 } else { /* incomplete Cholesky factor */ 245 factinfo.fill = 5.0; 246 factinfo.levels = lvl; 247 248 ierr = MatGetFactor(sA,MATSOLVERPETSC,MAT_FACTOR_ICC,&sC);CHKERRQ(ierr); 249 ierr = MatICCFactorSymbolic(sC,sA,perm,&factinfo);CHKERRQ(ierr); 250 } 251 ierr = MatCholeskyFactorNumeric(sC,sA,&factinfo);CHKERRQ(ierr); 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 ierr = MatDuplicate(sA,MAT_COPY_VALUES,&B);CHKERRQ(ierr); 257 ierr = MatICCFactor(B,perm,&factinfo);CHKERRQ(ierr); 258 ierr = MatEqual(sC,B,&equal);CHKERRQ(ierr); 259 if (!equal) { 260 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_USER,"in-place Cholesky factor != out-place Cholesky factor"); 261 } 262 ierr = MatDestroy(&B);CHKERRQ(ierr); 263 */ 264 } 265 266 ierr = MatMult(sA,x,b);CHKERRQ(ierr); 267 ierr = MatSolve(sC,b,y);CHKERRQ(ierr); 268 269 /* Test MatSolves() */ 270 if (bs == 1) { 271 Vecs xx,bb; 272 ierr = VecsCreateSeq(PETSC_COMM_SELF,n,4,&xx);CHKERRQ(ierr); 273 ierr = VecsDuplicate(xx,&bb);CHKERRQ(ierr); 274 ierr = MatSolves(sC,bb,xx);CHKERRQ(ierr); 275 ierr = VecsDestroy(xx);CHKERRQ(ierr); 276 ierr = VecsDestroy(bb);CHKERRQ(ierr); 277 } 278 ierr = MatDestroy(&sC);CHKERRQ(ierr); 279 280 /* Check the residual */ 281 ierr = VecAXPY(y,neg_one,x);CHKERRQ(ierr); 282 ierr = VecNorm(y,NORM_2,&norm2);CHKERRQ(ierr); 283 if (displ) { 284 ierr = PetscPrintf(PETSC_COMM_WORLD," lvl: %D, residual: %g\n", lvl,(double)norm2);CHKERRQ(ierr); 285 } 286 } 287 288 ierr = ISDestroy(&perm);CHKERRQ(ierr); 289 ierr = MatDestroy(&A);CHKERRQ(ierr); 290 ierr = MatDestroy(&sA);CHKERRQ(ierr); 291 ierr = VecDestroy(&x);CHKERRQ(ierr); 292 ierr = VecDestroy(&y);CHKERRQ(ierr); 293 ierr = VecDestroy(&b);CHKERRQ(ierr); 294 ierr = PetscRandomDestroy(&rdm);CHKERRQ(ierr); 295 296 ierr = PetscFinalize(); 297 return ierr; 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