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