1 2 static char help[] = "Tests the various sequential routines in MATSEQSBAIJ format.\n"; 3 4 #include <petscmat.h> 5 6 int main(int argc,char **args) 7 { 8 PetscMPIInt size; 9 Vec x,y,b,s1,s2; 10 Mat A; /* linear system matrix */ 11 Mat sA,sB,sFactor,B,C; /* symmetric matrices */ 12 PetscInt n,mbs=16,bs=1,nz=3,prob=1,i,j,k1,k2,col[3],lf,block, row,Ii,J,n1,inc; 13 PetscReal norm1,norm2,rnorm,tol=10*PETSC_SMALL; 14 PetscScalar neg_one=-1.0,four=4.0,value[3]; 15 IS perm, iscol; 16 PetscRandom rdm; 17 PetscBool doIcc=PETSC_TRUE,equal; 18 MatInfo minfo1,minfo2; 19 MatFactorInfo factinfo; 20 MatType type; 21 22 PetscCall(PetscInitialize(&argc,&args,(char*)0,help)); 23 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 24 PetscCheck(size == 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only!"); 25 PetscCall(PetscOptionsGetInt(NULL,NULL,"-bs",&bs,NULL)); 26 PetscCall(PetscOptionsGetInt(NULL,NULL,"-mbs",&mbs,NULL)); 27 28 n = mbs*bs; 29 PetscCall(MatCreate(PETSC_COMM_SELF,&A)); 30 PetscCall(MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE)); 31 PetscCall(MatSetType(A,MATSEQBAIJ)); 32 PetscCall(MatSetFromOptions(A)); 33 PetscCall(MatSeqBAIJSetPreallocation(A,bs,nz,NULL)); 34 35 PetscCall(MatCreate(PETSC_COMM_SELF,&sA)); 36 PetscCall(MatSetSizes(sA,n,n,PETSC_DETERMINE,PETSC_DETERMINE)); 37 PetscCall(MatSetType(sA,MATSEQSBAIJ)); 38 PetscCall(MatSetFromOptions(sA)); 39 PetscCall(MatGetType(sA,&type)); 40 PetscCall(PetscObjectTypeCompare((PetscObject)sA,MATSEQSBAIJ,&doIcc)); 41 PetscCall(MatSeqSBAIJSetPreallocation(sA,bs,nz,NULL)); 42 PetscCall(MatSetOption(sA,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_TRUE)); 43 44 /* Test MatGetOwnershipRange() */ 45 PetscCall(MatGetOwnershipRange(A,&Ii,&J)); 46 PetscCall(MatGetOwnershipRange(sA,&i,&j)); 47 if (i-Ii || j-J) { 48 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Error: MatGetOwnershipRange() in MatSBAIJ format\n")); 49 } 50 51 /* Assemble matrix */ 52 if (bs == 1) { 53 PetscCall(PetscOptionsGetInt(NULL,NULL,"-test_problem",&prob,NULL)); 54 if (prob == 1) { /* tridiagonal matrix */ 55 value[0] = -1.0; value[1] = 2.0; value[2] = -1.0; 56 for (i=1; i<n-1; i++) { 57 col[0] = i-1; col[1] = i; col[2] = i+1; 58 PetscCall(MatSetValues(A,1,&i,3,col,value,INSERT_VALUES)); 59 PetscCall(MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES)); 60 } 61 i = n - 1; col[0]=0; col[1] = n - 2; col[2] = n - 1; 62 63 value[0]= 0.1; value[1]=-1; value[2]=2; 64 65 PetscCall(MatSetValues(A,1,&i,3,col,value,INSERT_VALUES)); 66 PetscCall(MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES)); 67 68 i = 0; 69 col[0] = n-1; col[1] = 1; col[2] = 0; 70 value[0] = 0.1; value[1] = -1.0; value[2] = 2; 71 72 PetscCall(MatSetValues(A,1,&i,3,col,value,INSERT_VALUES)); 73 PetscCall(MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES)); 74 75 } else if (prob == 2) { /* matrix for the five point stencil */ 76 n1 = (PetscInt) (PetscSqrtReal((PetscReal)n) + 0.001); 77 PetscCheckFalse(n1*n1 - n,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"sqrt(n) must be a positive integer!"); 78 for (i=0; i<n1; i++) { 79 for (j=0; j<n1; j++) { 80 Ii = j + n1*i; 81 if (i>0) { 82 J = Ii - n1; 83 PetscCall(MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES)); 84 PetscCall(MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES)); 85 } 86 if (i<n1-1) { 87 J = Ii + n1; 88 PetscCall(MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES)); 89 PetscCall(MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES)); 90 } 91 if (j>0) { 92 J = Ii - 1; 93 PetscCall(MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES)); 94 PetscCall(MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES)); 95 } 96 if (j<n1-1) { 97 J = Ii + 1; 98 PetscCall(MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES)); 99 PetscCall(MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES)); 100 } 101 PetscCall(MatSetValues(A,1,&Ii,1,&Ii,&four,INSERT_VALUES)); 102 PetscCall(MatSetValues(sA,1,&Ii,1,&Ii,&four,INSERT_VALUES)); 103 } 104 } 105 } 106 107 } else { /* bs > 1 */ 108 for (block=0; block<n/bs; block++) { 109 /* diagonal blocks */ 110 value[0] = -1.0; value[1] = 4.0; value[2] = -1.0; 111 for (i=1+block*bs; i<bs-1+block*bs; i++) { 112 col[0] = i-1; col[1] = i; col[2] = i+1; 113 PetscCall(MatSetValues(A,1,&i,3,col,value,INSERT_VALUES)); 114 PetscCall(MatSetValues(sA,1,&i,3,col,value,INSERT_VALUES)); 115 } 116 i = bs - 1+block*bs; col[0] = bs - 2+block*bs; col[1] = bs - 1+block*bs; 117 118 value[0]=-1.0; value[1]=4.0; 119 120 PetscCall(MatSetValues(A,1,&i,2,col,value,INSERT_VALUES)); 121 PetscCall(MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES)); 122 123 i = 0+block*bs; col[0] = 0+block*bs; col[1] = 1+block*bs; 124 125 value[0]=4.0; value[1] = -1.0; 126 127 PetscCall(MatSetValues(A,1,&i,2,col,value,INSERT_VALUES)); 128 PetscCall(MatSetValues(sA,1,&i,2,col,value,INSERT_VALUES)); 129 } 130 /* off-diagonal blocks */ 131 value[0]=-1.0; 132 for (i=0; i<(n/bs-1)*bs; i++) { 133 col[0]=i+bs; 134 135 PetscCall(MatSetValues(A,1,&i,1,col,value,INSERT_VALUES)); 136 PetscCall(MatSetValues(sA,1,&i,1,col,value,INSERT_VALUES)); 137 138 col[0]=i; row=i+bs; 139 140 PetscCall(MatSetValues(A,1,&row,1,col,value,INSERT_VALUES)); 141 PetscCall(MatSetValues(sA,1,&row,1,col,value,INSERT_VALUES)); 142 } 143 } 144 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 145 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 146 147 PetscCall(MatAssemblyBegin(sA,MAT_FINAL_ASSEMBLY)); 148 PetscCall(MatAssemblyEnd(sA,MAT_FINAL_ASSEMBLY)); 149 150 /* Test MatGetInfo() of A and sA */ 151 PetscCall(MatGetInfo(A,MAT_LOCAL,&minfo1)); 152 PetscCall(MatGetInfo(sA,MAT_LOCAL,&minfo2)); 153 i = (int) (minfo1.nz_used - minfo2.nz_used); 154 j = (int) (minfo1.nz_allocated - minfo2.nz_allocated); 155 k1 = (int) (minfo1.nz_allocated - minfo1.nz_used); 156 k2 = (int) (minfo2.nz_allocated - minfo2.nz_used); 157 if (i < 0 || j < 0 || k1 < 0 || k2 < 0) { 158 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Error (compare A and sA): MatGetInfo()\n")); 159 } 160 161 /* Test MatDuplicate() */ 162 PetscCall(MatNorm(A,NORM_FROBENIUS,&norm1)); 163 PetscCall(MatDuplicate(sA,MAT_COPY_VALUES,&sB)); 164 PetscCall(MatEqual(sA,sB,&equal)); 165 PetscCheck(equal,PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"Error in MatDuplicate()"); 166 167 /* Test MatNorm() */ 168 PetscCall(MatNorm(A,NORM_FROBENIUS,&norm1)); 169 PetscCall(MatNorm(sB,NORM_FROBENIUS,&norm2)); 170 rnorm = PetscAbsReal(norm1-norm2)/norm2; 171 if (rnorm > tol) { 172 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_FROBENIUS, NormA=%16.14e NormsB=%16.14e\n",(double)norm1,(double)norm2)); 173 } 174 PetscCall(MatNorm(A,NORM_INFINITY,&norm1)); 175 PetscCall(MatNorm(sB,NORM_INFINITY,&norm2)); 176 rnorm = PetscAbsReal(norm1-norm2)/norm2; 177 if (rnorm > tol) { 178 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_INFINITY(), NormA=%16.14e NormsB=%16.14e\n",(double)norm1,(double)norm2)); 179 } 180 PetscCall(MatNorm(A,NORM_1,&norm1)); 181 PetscCall(MatNorm(sB,NORM_1,&norm2)); 182 rnorm = PetscAbsReal(norm1-norm2)/norm2; 183 if (rnorm > tol) { 184 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_INFINITY(), NormA=%16.14e NormsB=%16.14e\n",(double)norm1,(double)norm2)); 185 } 186 187 /* Test MatGetInfo(), MatGetSize(), MatGetBlockSize() */ 188 PetscCall(MatGetInfo(A,MAT_LOCAL,&minfo1)); 189 PetscCall(MatGetInfo(sB,MAT_LOCAL,&minfo2)); 190 i = (int) (minfo1.nz_used - minfo2.nz_used); 191 j = (int) (minfo1.nz_allocated - minfo2.nz_allocated); 192 k1 = (int) (minfo1.nz_allocated - minfo1.nz_used); 193 k2 = (int) (minfo2.nz_allocated - minfo2.nz_used); 194 if (i < 0 || j < 0 || k1 < 0 || k2 < 0) { 195 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Error(compare A and sB): MatGetInfo()\n")); 196 } 197 198 PetscCall(MatGetSize(A,&Ii,&J)); 199 PetscCall(MatGetSize(sB,&i,&j)); 200 if (i-Ii || j-J) { 201 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Error: MatGetSize()\n")); 202 } 203 204 PetscCall(MatGetBlockSize(A, &Ii)); 205 PetscCall(MatGetBlockSize(sB, &i)); 206 if (i-Ii) { 207 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Error: MatGetBlockSize()\n")); 208 } 209 210 PetscCall(PetscRandomCreate(PETSC_COMM_SELF,&rdm)); 211 PetscCall(PetscRandomSetFromOptions(rdm)); 212 PetscCall(VecCreateSeq(PETSC_COMM_SELF,n,&x)); 213 PetscCall(VecDuplicate(x,&s1)); 214 PetscCall(VecDuplicate(x,&s2)); 215 PetscCall(VecDuplicate(x,&y)); 216 PetscCall(VecDuplicate(x,&b)); 217 PetscCall(VecSetRandom(x,rdm)); 218 219 /* Test MatDiagonalScale(), MatGetDiagonal(), MatScale() */ 220 #if !defined(PETSC_USE_COMPLEX) 221 /* Scaling matrix with complex numbers results non-spd matrix, 222 causing crash of MatForwardSolve() and MatBackwardSolve() */ 223 PetscCall(MatDiagonalScale(A,x,x)); 224 PetscCall(MatDiagonalScale(sB,x,x)); 225 PetscCall(MatMultEqual(A,sB,10,&equal)); 226 PetscCheck(equal,PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"Error in MatDiagonalScale"); 227 228 PetscCall(MatGetDiagonal(A,s1)); 229 PetscCall(MatGetDiagonal(sB,s2)); 230 PetscCall(VecAXPY(s2,neg_one,s1)); 231 PetscCall(VecNorm(s2,NORM_1,&norm1)); 232 if (norm1>tol) { 233 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Error:MatGetDiagonal(), ||s1-s2||=%g\n",(double)norm1)); 234 } 235 236 { 237 PetscScalar alpha=0.1; 238 PetscCall(MatScale(A,alpha)); 239 PetscCall(MatScale(sB,alpha)); 240 } 241 #endif 242 243 /* Test MatGetRowMaxAbs() */ 244 PetscCall(MatGetRowMaxAbs(A,s1,NULL)); 245 PetscCall(MatGetRowMaxAbs(sB,s2,NULL)); 246 PetscCall(VecNorm(s1,NORM_1,&norm1)); 247 PetscCall(VecNorm(s2,NORM_1,&norm2)); 248 norm1 -= norm2; 249 if (norm1<-tol || norm1>tol) { 250 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Error:MatGetRowMaxAbs() \n")); 251 } 252 253 /* Test MatMult() */ 254 for (i=0; i<40; i++) { 255 PetscCall(VecSetRandom(x,rdm)); 256 PetscCall(MatMult(A,x,s1)); 257 PetscCall(MatMult(sB,x,s2)); 258 PetscCall(VecNorm(s1,NORM_1,&norm1)); 259 PetscCall(VecNorm(s2,NORM_1,&norm2)); 260 norm1 -= norm2; 261 if (norm1<-tol || norm1>tol) { 262 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Error: MatMult(), norm1-norm2: %g\n",(double)norm1)); 263 } 264 } 265 266 /* MatMultAdd() */ 267 for (i=0; i<40; i++) { 268 PetscCall(VecSetRandom(x,rdm)); 269 PetscCall(VecSetRandom(y,rdm)); 270 PetscCall(MatMultAdd(A,x,y,s1)); 271 PetscCall(MatMultAdd(sB,x,y,s2)); 272 PetscCall(VecNorm(s1,NORM_1,&norm1)); 273 PetscCall(VecNorm(s2,NORM_1,&norm2)); 274 norm1 -= norm2; 275 if (norm1<-tol || norm1>tol) { 276 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Error:MatMultAdd(), norm1-norm2: %g\n",(double)norm1)); 277 } 278 } 279 280 /* Test MatMatMult() for sbaij and dense matrices */ 281 PetscCall(MatCreateSeqDense(PETSC_COMM_SELF,n,5*n,NULL,&B)); 282 PetscCall(MatSetRandom(B,rdm)); 283 PetscCall(MatMatMult(sA,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C)); 284 PetscCall(MatMatMultEqual(sA,B,C,5*n,&equal)); 285 PetscCheck(equal,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error: MatMatMult()"); 286 PetscCall(MatDestroy(&C)); 287 PetscCall(MatDestroy(&B)); 288 289 /* Test MatCholeskyFactor(), MatICCFactor() with natural ordering */ 290 PetscCall(MatGetOrdering(A,MATORDERINGNATURAL,&perm,&iscol)); 291 PetscCall(ISDestroy(&iscol)); 292 norm1 = tol; 293 inc = bs; 294 295 /* initialize factinfo */ 296 PetscCall(PetscMemzero(&factinfo,sizeof(MatFactorInfo))); 297 298 for (lf=-1; lf<10; lf += inc) { 299 if (lf==-1) { /* Cholesky factor of sB (duplicate sA) */ 300 factinfo.fill = 5.0; 301 302 PetscCall(MatGetFactor(sB,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&sFactor)); 303 PetscCall(MatCholeskyFactorSymbolic(sFactor,sB,perm,&factinfo)); 304 } else if (!doIcc) break; 305 else { /* incomplete Cholesky factor */ 306 factinfo.fill = 5.0; 307 factinfo.levels = lf; 308 309 PetscCall(MatGetFactor(sB,MATSOLVERPETSC,MAT_FACTOR_ICC,&sFactor)); 310 PetscCall(MatICCFactorSymbolic(sFactor,sB,perm,&factinfo)); 311 } 312 PetscCall(MatCholeskyFactorNumeric(sFactor,sB,&factinfo)); 313 /* MatView(sFactor, PETSC_VIEWER_DRAW_WORLD); */ 314 315 /* test MatGetDiagonal on numeric factor */ 316 /* 317 if (lf == -1) { 318 PetscCall(MatGetDiagonal(sFactor,s1)); 319 printf(" in ex74.c, diag: \n"); 320 PetscCall(VecView(s1,PETSC_VIEWER_STDOUT_SELF)); 321 } 322 */ 323 324 PetscCall(MatMult(sB,x,b)); 325 326 /* test MatForwardSolve() and MatBackwardSolve() */ 327 if (lf == -1) { 328 PetscCall(MatForwardSolve(sFactor,b,s1)); 329 PetscCall(MatBackwardSolve(sFactor,s1,s2)); 330 PetscCall(VecAXPY(s2,neg_one,x)); 331 PetscCall(VecNorm(s2,NORM_2,&norm2)); 332 if (10*norm1 < norm2) { 333 PetscCall(PetscPrintf(PETSC_COMM_SELF,"MatForwardSolve and BackwardSolve: Norm of error=%g, bs=%" PetscInt_FMT "\n",(double)norm2,bs)); 334 } 335 } 336 337 /* test MatSolve() */ 338 PetscCall(MatSolve(sFactor,b,y)); 339 PetscCall(MatDestroy(&sFactor)); 340 /* Check the error */ 341 PetscCall(VecAXPY(y,neg_one,x)); 342 PetscCall(VecNorm(y,NORM_2,&norm2)); 343 if (10*norm1 < norm2 && lf-inc != -1) { 344 PetscCall(PetscPrintf(PETSC_COMM_SELF,"lf=%" PetscInt_FMT ", %" PetscInt_FMT ", Norm of error=%g, %g\n",lf-inc,lf,(double)norm1,(double)norm2)); 345 } 346 norm1 = norm2; 347 if (norm2 < tol && lf != -1) break; 348 } 349 350 #if defined(PETSC_HAVE_MUMPS) 351 PetscCall(MatGetFactor(sA,MATSOLVERMUMPS,MAT_FACTOR_CHOLESKY,&sFactor)); 352 PetscCall(MatCholeskyFactorSymbolic(sFactor,sA,NULL,NULL)); 353 PetscCall(MatCholeskyFactorNumeric(sFactor,sA,NULL)); 354 for (i=0; i<10; i++) { 355 PetscCall(VecSetRandom(b,rdm)); 356 PetscCall(MatSolve(sFactor,b,y)); 357 /* Check the error */ 358 PetscCall(MatMult(sA,y,x)); 359 PetscCall(VecAXPY(x,neg_one,b)); 360 PetscCall(VecNorm(x,NORM_2,&norm2)); 361 if (norm2>tol) { 362 PetscCall(PetscPrintf(PETSC_COMM_SELF,"Error:MatSolve(), norm2: %g\n",(double)norm2)); 363 } 364 } 365 PetscCall(MatDestroy(&sFactor)); 366 #endif 367 368 PetscCall(ISDestroy(&perm)); 369 370 PetscCall(MatDestroy(&A)); 371 PetscCall(MatDestroy(&sB)); 372 PetscCall(MatDestroy(&sA)); 373 PetscCall(VecDestroy(&x)); 374 PetscCall(VecDestroy(&y)); 375 PetscCall(VecDestroy(&s1)); 376 PetscCall(VecDestroy(&s2)); 377 PetscCall(VecDestroy(&b)); 378 PetscCall(PetscRandomDestroy(&rdm)); 379 380 PetscCall(PetscFinalize()); 381 return 0; 382 } 383 384 /*TEST 385 386 test: 387 args: -bs {{1 2 3 4 5 6 7 8}} 388 389 TEST*/ 390