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