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 PetscErrorCode ierr; 10 Vec x,y,b,s1,s2; 11 Mat A; /* linear system matrix */ 12 Mat sA,sB,sFactor,B,C; /* symmetric matrices */ 13 PetscInt n,mbs=16,bs=1,nz=3,prob=1,i,j,k1,k2,col[3],lf,block, row,Ii,J,n1,inc; 14 PetscReal norm1,norm2,rnorm,tol=10*PETSC_SMALL; 15 PetscScalar neg_one=-1.0,four=4.0,value[3]; 16 IS perm, iscol; 17 PetscRandom rdm; 18 PetscBool doIcc=PETSC_TRUE,equal; 19 MatInfo minfo1,minfo2; 20 MatFactorInfo factinfo; 21 MatType type; 22 23 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; 24 CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 25 PetscCheckFalse(size != 1,PETSC_COMM_WORLD,PETSC_ERR_SUP,"This is a uniprocessor example only!"); 26 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-bs",&bs,NULL)); 27 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-mbs",&mbs,NULL)); 28 29 n = mbs*bs; 30 CHKERRQ(MatCreate(PETSC_COMM_SELF,&A)); 31 CHKERRQ(MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE)); 32 CHKERRQ(MatSetType(A,MATSEQBAIJ)); 33 CHKERRQ(MatSetFromOptions(A)); 34 CHKERRQ(MatSeqBAIJSetPreallocation(A,bs,nz,NULL)); 35 36 CHKERRQ(MatCreate(PETSC_COMM_SELF,&sA)); 37 CHKERRQ(MatSetSizes(sA,n,n,PETSC_DETERMINE,PETSC_DETERMINE)); 38 CHKERRQ(MatSetType(sA,MATSEQSBAIJ)); 39 CHKERRQ(MatSetFromOptions(sA)); 40 CHKERRQ(MatGetType(sA,&type)); 41 CHKERRQ(PetscObjectTypeCompare((PetscObject)sA,MATSEQSBAIJ,&doIcc)); 42 CHKERRQ(MatSeqSBAIJSetPreallocation(sA,bs,nz,NULL)); 43 CHKERRQ(MatSetOption(sA,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_TRUE)); 44 45 /* Test MatGetOwnershipRange() */ 46 CHKERRQ(MatGetOwnershipRange(A,&Ii,&J)); 47 CHKERRQ(MatGetOwnershipRange(sA,&i,&j)); 48 if (i-Ii || j-J) { 49 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"Error: MatGetOwnershipRange() in MatSBAIJ format\n")); 50 } 51 52 /* Assemble matrix */ 53 if (bs == 1) { 54 CHKERRQ(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 CHKERRQ(MatSetValues(A,1,&i,3,col,value,INSERT_VALUES)); 60 CHKERRQ(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 CHKERRQ(MatSetValues(A,1,&i,3,col,value,INSERT_VALUES)); 67 CHKERRQ(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 CHKERRQ(MatSetValues(A,1,&i,3,col,value,INSERT_VALUES)); 74 CHKERRQ(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 PetscCheckFalse(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 CHKERRQ(MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES)); 85 CHKERRQ(MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES)); 86 } 87 if (i<n1-1) { 88 J = Ii + n1; 89 CHKERRQ(MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES)); 90 CHKERRQ(MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES)); 91 } 92 if (j>0) { 93 J = Ii - 1; 94 CHKERRQ(MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES)); 95 CHKERRQ(MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES)); 96 } 97 if (j<n1-1) { 98 J = Ii + 1; 99 CHKERRQ(MatSetValues(A,1,&Ii,1,&J,&neg_one,INSERT_VALUES)); 100 CHKERRQ(MatSetValues(sA,1,&Ii,1,&J,&neg_one,INSERT_VALUES)); 101 } 102 CHKERRQ(MatSetValues(A,1,&Ii,1,&Ii,&four,INSERT_VALUES)); 103 CHKERRQ(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 CHKERRQ(MatSetValues(A,1,&i,3,col,value,INSERT_VALUES)); 115 CHKERRQ(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 CHKERRQ(MatSetValues(A,1,&i,2,col,value,INSERT_VALUES)); 122 CHKERRQ(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 CHKERRQ(MatSetValues(A,1,&i,2,col,value,INSERT_VALUES)); 129 CHKERRQ(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 CHKERRQ(MatSetValues(A,1,&i,1,col,value,INSERT_VALUES)); 137 CHKERRQ(MatSetValues(sA,1,&i,1,col,value,INSERT_VALUES)); 138 139 col[0]=i; row=i+bs; 140 141 CHKERRQ(MatSetValues(A,1,&row,1,col,value,INSERT_VALUES)); 142 CHKERRQ(MatSetValues(sA,1,&row,1,col,value,INSERT_VALUES)); 143 } 144 } 145 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 146 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 147 148 CHKERRQ(MatAssemblyBegin(sA,MAT_FINAL_ASSEMBLY)); 149 CHKERRQ(MatAssemblyEnd(sA,MAT_FINAL_ASSEMBLY)); 150 151 /* Test MatGetInfo() of A and sA */ 152 CHKERRQ(MatGetInfo(A,MAT_LOCAL,&minfo1)); 153 CHKERRQ(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 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"Error (compare A and sA): MatGetInfo()\n")); 160 } 161 162 /* Test MatDuplicate() */ 163 CHKERRQ(MatNorm(A,NORM_FROBENIUS,&norm1)); 164 CHKERRQ(MatDuplicate(sA,MAT_COPY_VALUES,&sB)); 165 CHKERRQ(MatEqual(sA,sB,&equal)); 166 PetscCheck(equal,PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"Error in MatDuplicate()"); 167 168 /* Test MatNorm() */ 169 CHKERRQ(MatNorm(A,NORM_FROBENIUS,&norm1)); 170 CHKERRQ(MatNorm(sB,NORM_FROBENIUS,&norm2)); 171 rnorm = PetscAbsReal(norm1-norm2)/norm2; 172 if (rnorm > tol) { 173 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_FROBENIUS, NormA=%16.14e NormsB=%16.14e\n",(double)norm1,(double)norm2)); 174 } 175 CHKERRQ(MatNorm(A,NORM_INFINITY,&norm1)); 176 CHKERRQ(MatNorm(sB,NORM_INFINITY,&norm2)); 177 rnorm = PetscAbsReal(norm1-norm2)/norm2; 178 if (rnorm > tol) { 179 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"Error: MatNorm_INFINITY(), NormA=%16.14e NormsB=%16.14e\n",(double)norm1,(double)norm2)); 180 } 181 CHKERRQ(MatNorm(A,NORM_1,&norm1)); 182 CHKERRQ(MatNorm(sB,NORM_1,&norm2)); 183 rnorm = PetscAbsReal(norm1-norm2)/norm2; 184 if (rnorm > tol) { 185 CHKERRQ(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 CHKERRQ(MatGetInfo(A,MAT_LOCAL,&minfo1)); 190 CHKERRQ(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 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"Error(compare A and sB): MatGetInfo()\n")); 197 } 198 199 CHKERRQ(MatGetSize(A,&Ii,&J)); 200 CHKERRQ(MatGetSize(sB,&i,&j)); 201 if (i-Ii || j-J) { 202 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"Error: MatGetSize()\n")); 203 } 204 205 CHKERRQ(MatGetBlockSize(A, &Ii)); 206 CHKERRQ(MatGetBlockSize(sB, &i)); 207 if (i-Ii) { 208 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"Error: MatGetBlockSize()\n")); 209 } 210 211 CHKERRQ(PetscRandomCreate(PETSC_COMM_SELF,&rdm)); 212 CHKERRQ(PetscRandomSetFromOptions(rdm)); 213 CHKERRQ(VecCreateSeq(PETSC_COMM_SELF,n,&x)); 214 CHKERRQ(VecDuplicate(x,&s1)); 215 CHKERRQ(VecDuplicate(x,&s2)); 216 CHKERRQ(VecDuplicate(x,&y)); 217 CHKERRQ(VecDuplicate(x,&b)); 218 CHKERRQ(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 CHKERRQ(MatDiagonalScale(A,x,x)); 225 CHKERRQ(MatDiagonalScale(sB,x,x)); 226 CHKERRQ(MatMultEqual(A,sB,10,&equal)); 227 PetscCheck(equal,PETSC_COMM_SELF,PETSC_ERR_ARG_NOTSAMETYPE,"Error in MatDiagonalScale"); 228 229 CHKERRQ(MatGetDiagonal(A,s1)); 230 CHKERRQ(MatGetDiagonal(sB,s2)); 231 CHKERRQ(VecAXPY(s2,neg_one,s1)); 232 CHKERRQ(VecNorm(s2,NORM_1,&norm1)); 233 if (norm1>tol) { 234 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"Error:MatGetDiagonal(), ||s1-s2||=%g\n",(double)norm1)); 235 } 236 237 { 238 PetscScalar alpha=0.1; 239 CHKERRQ(MatScale(A,alpha)); 240 CHKERRQ(MatScale(sB,alpha)); 241 } 242 #endif 243 244 /* Test MatGetRowMaxAbs() */ 245 CHKERRQ(MatGetRowMaxAbs(A,s1,NULL)); 246 CHKERRQ(MatGetRowMaxAbs(sB,s2,NULL)); 247 CHKERRQ(VecNorm(s1,NORM_1,&norm1)); 248 CHKERRQ(VecNorm(s2,NORM_1,&norm2)); 249 norm1 -= norm2; 250 if (norm1<-tol || norm1>tol) { 251 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"Error:MatGetRowMaxAbs() \n")); 252 } 253 254 /* Test MatMult() */ 255 for (i=0; i<40; i++) { 256 CHKERRQ(VecSetRandom(x,rdm)); 257 CHKERRQ(MatMult(A,x,s1)); 258 CHKERRQ(MatMult(sB,x,s2)); 259 CHKERRQ(VecNorm(s1,NORM_1,&norm1)); 260 CHKERRQ(VecNorm(s2,NORM_1,&norm2)); 261 norm1 -= norm2; 262 if (norm1<-tol || norm1>tol) { 263 CHKERRQ(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 CHKERRQ(VecSetRandom(x,rdm)); 270 CHKERRQ(VecSetRandom(y,rdm)); 271 CHKERRQ(MatMultAdd(A,x,y,s1)); 272 CHKERRQ(MatMultAdd(sB,x,y,s2)); 273 CHKERRQ(VecNorm(s1,NORM_1,&norm1)); 274 CHKERRQ(VecNorm(s2,NORM_1,&norm2)); 275 norm1 -= norm2; 276 if (norm1<-tol || norm1>tol) { 277 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"Error:MatMultAdd(), norm1-norm2: %g\n",(double)norm1)); 278 } 279 } 280 281 /* Test MatMatMult() for sbaij and dense matrices */ 282 CHKERRQ(MatCreateSeqDense(PETSC_COMM_SELF,n,5*n,NULL,&B)); 283 CHKERRQ(MatSetRandom(B,rdm)); 284 CHKERRQ(MatMatMult(sA,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C)); 285 CHKERRQ(MatMatMultEqual(sA,B,C,5*n,&equal)); 286 PetscCheck(equal,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error: MatMatMult()"); 287 CHKERRQ(MatDestroy(&C)); 288 CHKERRQ(MatDestroy(&B)); 289 290 /* Test MatCholeskyFactor(), MatICCFactor() with natural ordering */ 291 CHKERRQ(MatGetOrdering(A,MATORDERINGNATURAL,&perm,&iscol)); 292 CHKERRQ(ISDestroy(&iscol)); 293 norm1 = tol; 294 inc = bs; 295 296 /* initialize factinfo */ 297 CHKERRQ(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 CHKERRQ(MatGetFactor(sB,MATSOLVERPETSC,MAT_FACTOR_CHOLESKY,&sFactor)); 304 CHKERRQ(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 CHKERRQ(MatGetFactor(sB,MATSOLVERPETSC,MAT_FACTOR_ICC,&sFactor)); 311 CHKERRQ(MatICCFactorSymbolic(sFactor,sB,perm,&factinfo)); 312 } 313 CHKERRQ(MatCholeskyFactorNumeric(sFactor,sB,&factinfo)); 314 /* MatView(sFactor, PETSC_VIEWER_DRAW_WORLD); */ 315 316 /* test MatGetDiagonal on numeric factor */ 317 /* 318 if (lf == -1) { 319 CHKERRQ(MatGetDiagonal(sFactor,s1)); 320 printf(" in ex74.c, diag: \n"); 321 CHKERRQ(VecView(s1,PETSC_VIEWER_STDOUT_SELF)); 322 } 323 */ 324 325 CHKERRQ(MatMult(sB,x,b)); 326 327 /* test MatForwardSolve() and MatBackwardSolve() */ 328 if (lf == -1) { 329 CHKERRQ(MatForwardSolve(sFactor,b,s1)); 330 CHKERRQ(MatBackwardSolve(sFactor,s1,s2)); 331 CHKERRQ(VecAXPY(s2,neg_one,x)); 332 CHKERRQ(VecNorm(s2,NORM_2,&norm2)); 333 if (10*norm1 < norm2) { 334 CHKERRQ(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 CHKERRQ(MatSolve(sFactor,b,y)); 340 CHKERRQ(MatDestroy(&sFactor)); 341 /* Check the error */ 342 CHKERRQ(VecAXPY(y,neg_one,x)); 343 CHKERRQ(VecNorm(y,NORM_2,&norm2)); 344 if (10*norm1 < norm2 && lf-inc != -1) { 345 CHKERRQ(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 CHKERRQ(MatGetFactor(sA,MATSOLVERMUMPS,MAT_FACTOR_CHOLESKY,&sFactor)); 353 CHKERRQ(MatCholeskyFactorSymbolic(sFactor,sA,NULL,NULL)); 354 CHKERRQ(MatCholeskyFactorNumeric(sFactor,sA,NULL)); 355 for (i=0; i<10; i++) { 356 CHKERRQ(VecSetRandom(b,rdm)); 357 CHKERRQ(MatSolve(sFactor,b,y)); 358 /* Check the error */ 359 CHKERRQ(MatMult(sA,y,x)); 360 CHKERRQ(VecAXPY(x,neg_one,b)); 361 CHKERRQ(VecNorm(x,NORM_2,&norm2)); 362 if (norm2>tol) { 363 CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"Error:MatSolve(), norm2: %g\n",(double)norm2)); 364 } 365 } 366 CHKERRQ(MatDestroy(&sFactor)); 367 #endif 368 369 CHKERRQ(ISDestroy(&perm)); 370 371 CHKERRQ(MatDestroy(&A)); 372 CHKERRQ(MatDestroy(&sB)); 373 CHKERRQ(MatDestroy(&sA)); 374 CHKERRQ(VecDestroy(&x)); 375 CHKERRQ(VecDestroy(&y)); 376 CHKERRQ(VecDestroy(&s1)); 377 CHKERRQ(VecDestroy(&s2)); 378 CHKERRQ(VecDestroy(&b)); 379 CHKERRQ(PetscRandomDestroy(&rdm)); 380 381 ierr = PetscFinalize(); 382 return ierr; 383 } 384 385 /*TEST 386 387 test: 388 args: -bs {{1 2 3 4 5 6 7 8}} 389 390 TEST*/ 391