1 static char help[] = "Tests various routines for MATSHELL\n\n"; 2 3 #include <petscmat.h> 4 5 typedef struct _n_User *User; 6 struct _n_User { 7 Mat B; 8 }; 9 10 static PetscErrorCode MatGetDiagonal_User(Mat A,Vec X) 11 { 12 User user; 13 14 PetscFunctionBegin; 15 PetscCall(MatShellGetContext(A,&user)); 16 PetscCall(MatGetDiagonal(user->B,X)); 17 PetscFunctionReturn(0); 18 } 19 20 static PetscErrorCode MatMult_User(Mat A,Vec X,Vec Y) 21 { 22 User user; 23 24 PetscFunctionBegin; 25 PetscCall(MatShellGetContext(A,&user)); 26 PetscCall(MatMult(user->B,X,Y)); 27 PetscFunctionReturn(0); 28 } 29 30 static PetscErrorCode MatMultTranspose_User(Mat A,Vec X,Vec Y) 31 { 32 User user; 33 34 PetscFunctionBegin; 35 PetscCall(MatShellGetContext(A,&user)); 36 PetscCall(MatMultTranspose(user->B,X,Y)); 37 PetscFunctionReturn(0); 38 } 39 40 static PetscErrorCode MatCopy_User(Mat A,Mat X,MatStructure str) 41 { 42 User user,userX; 43 44 PetscFunctionBegin; 45 PetscCall(MatShellGetContext(A,&user)); 46 PetscCall(MatShellGetContext(X,&userX)); 47 PetscCheck(user == userX,PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"This should not happen"); 48 PetscCall(PetscObjectReference((PetscObject)user->B)); 49 PetscFunctionReturn(0); 50 } 51 52 static PetscErrorCode MatDestroy_User(Mat A) 53 { 54 User user; 55 56 PetscFunctionBegin; 57 PetscCall(MatShellGetContext(A,&user)); 58 PetscCall(PetscObjectDereference((PetscObject)user->B)); 59 PetscFunctionReturn(0); 60 } 61 62 int main(int argc,char **args) 63 { 64 User user; 65 Mat A,S; 66 PetscScalar *data,diag = 1.3; 67 PetscReal tol = PETSC_SMALL; 68 PetscInt i,j,m = PETSC_DECIDE,n = PETSC_DECIDE,M = 17,N = 15,s1,s2; 69 PetscInt test, ntest = 2; 70 PetscMPIInt rank,size; 71 PetscBool nc = PETSC_FALSE, cong, flg; 72 PetscBool ronl = PETSC_TRUE; 73 PetscBool randomize = PETSC_FALSE, submat = PETSC_FALSE; 74 PetscBool keep = PETSC_FALSE; 75 PetscBool testzerorows = PETSC_TRUE, testdiagscale = PETSC_TRUE, testgetdiag = PETSC_TRUE, testsubmat = PETSC_TRUE; 76 PetscBool testshift = PETSC_TRUE, testscale = PETSC_TRUE, testdup = PETSC_TRUE, testreset = PETSC_TRUE; 77 PetscBool testaxpy = PETSC_TRUE, testaxpyd = PETSC_TRUE, testaxpyerr = PETSC_FALSE; 78 79 PetscCall(PetscInitialize(&argc,&args,(char*)0,help)); 80 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 81 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 82 PetscCall(PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL)); 83 PetscCall(PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL)); 84 PetscCall(PetscOptionsGetInt(NULL,NULL,"-ml",&m,NULL)); 85 PetscCall(PetscOptionsGetInt(NULL,NULL,"-nl",&n,NULL)); 86 PetscCall(PetscOptionsGetBool(NULL,NULL,"-square_nc",&nc,NULL)); 87 PetscCall(PetscOptionsGetBool(NULL,NULL,"-rows_only",&ronl,NULL)); 88 PetscCall(PetscOptionsGetBool(NULL,NULL,"-randomize",&randomize,NULL)); 89 PetscCall(PetscOptionsGetBool(NULL,NULL,"-submat",&submat,NULL)); 90 PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_zerorows",&testzerorows,NULL)); 91 PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_diagscale",&testdiagscale,NULL)); 92 PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_getdiag",&testgetdiag,NULL)); 93 PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_shift",&testshift,NULL)); 94 PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_scale",&testscale,NULL)); 95 PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_dup",&testdup,NULL)); 96 PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_reset",&testreset,NULL)); 97 PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_submat",&testsubmat,NULL)); 98 PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_axpy",&testaxpy,NULL)); 99 PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_axpy_different",&testaxpyd,NULL)); 100 PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_axpy_error",&testaxpyerr,NULL)); 101 PetscCall(PetscOptionsGetInt(NULL,NULL,"-loop",&ntest,NULL)); 102 PetscCall(PetscOptionsGetReal(NULL,NULL,"-tol",&tol,NULL)); 103 PetscCall(PetscOptionsGetScalar(NULL,NULL,"-diag",&diag,NULL)); 104 PetscCall(PetscOptionsGetBool(NULL,NULL,"-keep",&keep,NULL)); 105 /* This tests square matrices with different row/col layout */ 106 if (nc && size > 1) { 107 M = PetscMax(PetscMax(N,M),1); 108 N = M; 109 m = n = 0; 110 if (rank == 0) { m = M-1; n = 1; } 111 else if (rank == 1) { m = 1; n = N-1; } 112 } 113 PetscCall(MatCreateDense(PETSC_COMM_WORLD,m,n,M,N,NULL,&A)); 114 PetscCall(MatGetLocalSize(A,&m,&n)); 115 PetscCall(MatGetSize(A,&M,&N)); 116 PetscCall(MatGetOwnershipRange(A,&s1,NULL)); 117 s2 = 1; 118 while (s2 < M) s2 *= 10; 119 PetscCall(MatDenseGetArray(A,&data)); 120 for (j = 0; j < N; j++) { 121 for (i = 0; i < m; i++) { 122 data[j*m + i] = s2*j + i + s1 + 1; 123 } 124 } 125 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 126 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 127 128 if (submat) { 129 Mat A2; 130 IS r,c; 131 PetscInt rst,ren,cst,cen; 132 133 PetscCall(MatGetOwnershipRange(A,&rst,&ren)); 134 PetscCall(MatGetOwnershipRangeColumn(A,&cst,&cen)); 135 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A),(ren-rst)/2,rst,1,&r)); 136 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A),(cen-cst)/2,cst,1,&c)); 137 PetscCall(MatCreateSubMatrix(A,r,c,MAT_INITIAL_MATRIX,&A2)); 138 PetscCall(ISDestroy(&r)); 139 PetscCall(ISDestroy(&c)); 140 PetscCall(MatDestroy(&A)); 141 A = A2; 142 } 143 144 PetscCall(MatGetSize(A,&M,&N)); 145 PetscCall(MatGetLocalSize(A,&m,&n)); 146 PetscCall(MatHasCongruentLayouts(A,&cong)); 147 148 PetscCall(MatConvert(A,MATAIJ,MAT_INPLACE_MATRIX,&A)); 149 PetscCall(MatSetOption(A,MAT_KEEP_NONZERO_PATTERN,keep)); 150 PetscCall(PetscObjectSetName((PetscObject)A,"initial")); 151 PetscCall(MatViewFromOptions(A,NULL,"-view_mat")); 152 153 PetscCall(PetscNew(&user)); 154 PetscCall(MatCreateShell(PETSC_COMM_WORLD,m,n,M,N,user,&S)); 155 PetscCall(MatShellSetOperation(S,MATOP_MULT,(void (*)(void))MatMult_User)); 156 PetscCall(MatShellSetOperation(S,MATOP_MULT_TRANSPOSE,(void (*)(void))MatMultTranspose_User)); 157 if (cong) PetscCall(MatShellSetOperation(S,MATOP_GET_DIAGONAL,(void (*)(void))MatGetDiagonal_User)); 158 PetscCall(MatShellSetOperation(S,MATOP_COPY,(void (*)(void))MatCopy_User)); 159 PetscCall(MatShellSetOperation(S,MATOP_DESTROY,(void (*)(void))MatDestroy_User)); 160 PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&user->B)); 161 162 /* Square and rows only scaling */ 163 ronl = cong ? ronl : PETSC_TRUE; 164 165 for (test = 0; test < ntest; test++) { 166 PetscReal err; 167 168 PetscCall(MatMultAddEqual(A,S,10,&flg)); 169 if (!flg) { 170 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error mult add\n",test)); 171 } 172 PetscCall(MatMultTransposeAddEqual(A,S,10,&flg)); 173 if (!flg) { 174 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error mult add (T)\n",test)); 175 } 176 if (testzerorows) { 177 Mat ST,B,C,BT,BTT; 178 IS zr; 179 Vec x = NULL, b1 = NULL, b2 = NULL; 180 PetscInt *idxs = NULL, nr = 0; 181 182 if (rank == (test%size)) { 183 nr = 1; 184 PetscCall(PetscMalloc1(nr,&idxs)); 185 if (test%2) { 186 idxs[0] = (2*M - 1 - test/2)%M; 187 } else { 188 idxs[0] = (test/2)%M; 189 } 190 idxs[0] = PetscMax(idxs[0],0); 191 } 192 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD,nr,idxs,PETSC_OWN_POINTER,&zr)); 193 PetscCall(PetscObjectSetName((PetscObject)zr,"ZR")); 194 PetscCall(ISViewFromOptions(zr,NULL,"-view_is")); 195 PetscCall(MatCreateVecs(A,&x,&b1)); 196 if (randomize) { 197 PetscCall(VecSetRandom(x,NULL)); 198 PetscCall(VecSetRandom(b1,NULL)); 199 } else { 200 PetscCall(VecSet(x,11.4)); 201 PetscCall(VecSet(b1,-14.2)); 202 } 203 PetscCall(VecDuplicate(b1,&b2)); 204 PetscCall(VecCopy(b1,b2)); 205 PetscCall(PetscObjectSetName((PetscObject)b1,"A_B1")); 206 PetscCall(PetscObjectSetName((PetscObject)b2,"A_B2")); 207 if (size > 1 && !cong) { /* MATMPIAIJ ZeroRows and ZeroRowsColumns are buggy in this case */ 208 PetscCall(VecDestroy(&b1)); 209 } 210 if (ronl) { 211 PetscCall(MatZeroRowsIS(A,zr,diag,x,b1)); 212 PetscCall(MatZeroRowsIS(S,zr,diag,x,b2)); 213 } else { 214 PetscCall(MatZeroRowsColumnsIS(A,zr,diag,x,b1)); 215 PetscCall(MatZeroRowsColumnsIS(S,zr,diag,x,b2)); 216 PetscCall(ISDestroy(&zr)); 217 /* Mix zerorows and zerorowscols */ 218 nr = 0; 219 idxs = NULL; 220 if (rank == 0) { 221 nr = 1; 222 PetscCall(PetscMalloc1(nr,&idxs)); 223 if (test%2) { 224 idxs[0] = (3*M - 2 - test/2)%M; 225 } else { 226 idxs[0] = (test/2+1)%M; 227 } 228 idxs[0] = PetscMax(idxs[0],0); 229 } 230 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD,nr,idxs,PETSC_OWN_POINTER,&zr)); 231 PetscCall(PetscObjectSetName((PetscObject)zr,"ZR2")); 232 PetscCall(ISViewFromOptions(zr,NULL,"-view_is")); 233 PetscCall(MatZeroRowsIS(A,zr,diag*2.0+PETSC_SMALL,NULL,NULL)); 234 PetscCall(MatZeroRowsIS(S,zr,diag*2.0+PETSC_SMALL,NULL,NULL)); 235 } 236 PetscCall(ISDestroy(&zr)); 237 238 if (b1) { 239 Vec b; 240 241 PetscCall(VecViewFromOptions(b1,NULL,"-view_b")); 242 PetscCall(VecViewFromOptions(b2,NULL,"-view_b")); 243 PetscCall(VecDuplicate(b1,&b)); 244 PetscCall(VecCopy(b1,b)); 245 PetscCall(VecAXPY(b,-1.0,b2)); 246 PetscCall(VecNorm(b,NORM_INFINITY,&err)); 247 if (err >= tol) { 248 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error b %g\n",test,(double)err)); 249 } 250 PetscCall(VecDestroy(&b)); 251 } 252 PetscCall(VecDestroy(&b1)); 253 PetscCall(VecDestroy(&b2)); 254 PetscCall(VecDestroy(&x)); 255 PetscCall(MatConvert(S,MATDENSE,MAT_INITIAL_MATRIX,&B)); 256 257 PetscCall(MatCreateTranspose(S,&ST)); 258 PetscCall(MatComputeOperator(ST,MATDENSE,&BT)); 259 PetscCall(MatTranspose(BT,MAT_INITIAL_MATRIX,&BTT)); 260 PetscCall(PetscObjectSetName((PetscObject)B,"S")); 261 PetscCall(PetscObjectSetName((PetscObject)BTT,"STT")); 262 PetscCall(MatConvert(A,MATDENSE,MAT_INITIAL_MATRIX,&C)); 263 PetscCall(PetscObjectSetName((PetscObject)C,"A")); 264 265 PetscCall(MatViewFromOptions(C,NULL,"-view_mat")); 266 PetscCall(MatViewFromOptions(B,NULL,"-view_mat")); 267 PetscCall(MatViewFromOptions(BTT,NULL,"-view_mat")); 268 269 PetscCall(MatAXPY(C,-1.0,B,SAME_NONZERO_PATTERN)); 270 PetscCall(MatNorm(C,NORM_FROBENIUS,&err)); 271 if (err >= tol) { 272 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error mat mult after %s %g\n",test,ronl ? "MatZeroRows" : "MatZeroRowsColumns",(double)err)); 273 } 274 275 PetscCall(MatConvert(A,MATDENSE,MAT_REUSE_MATRIX,&C)); 276 PetscCall(MatAXPY(C,-1.0,BTT,SAME_NONZERO_PATTERN)); 277 PetscCall(MatNorm(C,NORM_FROBENIUS,&err)); 278 if (err >= tol) { 279 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error mat mult transpose after %s %g\n",test,ronl ? "MatZeroRows" : "MatZeroRowsColumns",(double)err)); 280 } 281 282 PetscCall(MatDestroy(&ST)); 283 PetscCall(MatDestroy(&BTT)); 284 PetscCall(MatDestroy(&BT)); 285 PetscCall(MatDestroy(&B)); 286 PetscCall(MatDestroy(&C)); 287 } 288 if (testdiagscale) { /* MatDiagonalScale() */ 289 Vec vr,vl; 290 291 PetscCall(MatCreateVecs(A,&vr,&vl)); 292 if (randomize) { 293 PetscCall(VecSetRandom(vr,NULL)); 294 PetscCall(VecSetRandom(vl,NULL)); 295 } else { 296 PetscCall(VecSet(vr,test%2 ? 0.15 : 1.0/0.15)); 297 PetscCall(VecSet(vl,test%2 ? -1.2 : 1.0/-1.2)); 298 } 299 PetscCall(MatDiagonalScale(A,vl,vr)); 300 PetscCall(MatDiagonalScale(S,vl,vr)); 301 PetscCall(VecDestroy(&vr)); 302 PetscCall(VecDestroy(&vl)); 303 } 304 305 if (testscale) { /* MatScale() */ 306 PetscCall(MatScale(A,test%2 ? 1.4 : 1.0/1.4)); 307 PetscCall(MatScale(S,test%2 ? 1.4 : 1.0/1.4)); 308 } 309 310 if (testshift && cong) { /* MatShift() : MATSHELL shift is broken when row/cols layout are not congruent and left/right scaling have been applied */ 311 PetscCall(MatShift(A,test%2 ? -77.5 : 77.5)); 312 PetscCall(MatShift(S,test%2 ? -77.5 : 77.5)); 313 } 314 315 if (testgetdiag && cong) { /* MatGetDiagonal() */ 316 Vec dA,dS; 317 318 PetscCall(MatCreateVecs(A,&dA,NULL)); 319 PetscCall(MatCreateVecs(S,&dS,NULL)); 320 PetscCall(MatGetDiagonal(A,dA)); 321 PetscCall(MatGetDiagonal(S,dS)); 322 PetscCall(VecAXPY(dA,-1.0,dS)); 323 PetscCall(VecNorm(dA,NORM_INFINITY,&err)); 324 if (err >= tol) { 325 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error diag %g\n",test,(double)err)); 326 } 327 PetscCall(VecDestroy(&dA)); 328 PetscCall(VecDestroy(&dS)); 329 } 330 331 if (testdup && !test) { 332 Mat A2, S2; 333 334 PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&A2)); 335 PetscCall(MatDuplicate(S,MAT_COPY_VALUES,&S2)); 336 PetscCall(MatDestroy(&A)); 337 PetscCall(MatDestroy(&S)); 338 A = A2; 339 S = S2; 340 } 341 342 if (testsubmat) { 343 Mat sA,sS,dA,dS,At,St; 344 IS r,c; 345 PetscInt rst,ren,cst,cen; 346 347 PetscCall(MatGetOwnershipRange(A,&rst,&ren)); 348 PetscCall(MatGetOwnershipRangeColumn(A,&cst,&cen)); 349 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A),(ren-rst)/2,rst,1,&r)); 350 PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A),(cen-cst)/2,cst,1,&c)); 351 PetscCall(MatCreateSubMatrix(A,r,c,MAT_INITIAL_MATRIX,&sA)); 352 PetscCall(MatCreateSubMatrix(S,r,c,MAT_INITIAL_MATRIX,&sS)); 353 PetscCall(MatMultAddEqual(sA,sS,10,&flg)); 354 if (!flg) { 355 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error submatrix mult add\n",test)); 356 } 357 PetscCall(MatMultTransposeAddEqual(sA,sS,10,&flg)); 358 if (!flg) { 359 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error submatrix mult add (T)\n",test)); 360 } 361 PetscCall(MatConvert(sA,MATDENSE,MAT_INITIAL_MATRIX,&dA)); 362 PetscCall(MatConvert(sS,MATDENSE,MAT_INITIAL_MATRIX,&dS)); 363 PetscCall(MatAXPY(dA,-1.0,dS,SAME_NONZERO_PATTERN)); 364 PetscCall(MatNorm(dA,NORM_FROBENIUS,&err)); 365 if (err >= tol) { 366 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error mat submatrix %g\n",test,(double)err)); 367 } 368 PetscCall(MatDestroy(&sA)); 369 PetscCall(MatDestroy(&sS)); 370 PetscCall(MatDestroy(&dA)); 371 PetscCall(MatDestroy(&dS)); 372 PetscCall(MatCreateTranspose(A,&At)); 373 PetscCall(MatCreateTranspose(S,&St)); 374 PetscCall(MatCreateSubMatrix(At,c,r,MAT_INITIAL_MATRIX,&sA)); 375 PetscCall(MatCreateSubMatrix(St,c,r,MAT_INITIAL_MATRIX,&sS)); 376 PetscCall(MatMultAddEqual(sA,sS,10,&flg)); 377 if (!flg) { 378 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error submatrix (T) mult add\n",test)); 379 } 380 PetscCall(MatMultTransposeAddEqual(sA,sS,10,&flg)); 381 if (!flg) { 382 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error submatrix (T) mult add (T)\n",test)); 383 } 384 PetscCall(MatConvert(sA,MATDENSE,MAT_INITIAL_MATRIX,&dA)); 385 PetscCall(MatConvert(sS,MATDENSE,MAT_INITIAL_MATRIX,&dS)); 386 PetscCall(MatAXPY(dA,-1.0,dS,SAME_NONZERO_PATTERN)); 387 PetscCall(MatNorm(dA,NORM_FROBENIUS,&err)); 388 if (err >= tol) { 389 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error mat submatrix (T) %g\n",test,(double)err)); 390 } 391 PetscCall(MatDestroy(&sA)); 392 PetscCall(MatDestroy(&sS)); 393 PetscCall(MatDestroy(&dA)); 394 PetscCall(MatDestroy(&dS)); 395 PetscCall(MatDestroy(&At)); 396 PetscCall(MatDestroy(&St)); 397 PetscCall(ISDestroy(&r)); 398 PetscCall(ISDestroy(&c)); 399 } 400 401 if (testaxpy) { 402 Mat tA,tS,dA,dS; 403 MatStructure str[3] = { SAME_NONZERO_PATTERN, SUBSET_NONZERO_PATTERN, DIFFERENT_NONZERO_PATTERN }; 404 405 PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&tA)); 406 if (testaxpyd && !(test%2)) { 407 PetscCall(PetscObjectReference((PetscObject)tA)); 408 tS = tA; 409 } else { 410 PetscCall(PetscObjectReference((PetscObject)S)); 411 tS = S; 412 } 413 PetscCall(MatAXPY(A,0.5,tA,str[test%3])); 414 PetscCall(MatAXPY(S,0.5,tS,str[test%3])); 415 /* this will trigger an error the next MatMult or MatMultTranspose call for S */ 416 if (testaxpyerr) PetscCall(MatScale(tA,0)); 417 PetscCall(MatDestroy(&tA)); 418 PetscCall(MatDestroy(&tS)); 419 PetscCall(MatMultAddEqual(A,S,10,&flg)); 420 if (!flg) { 421 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error axpy mult add\n",test)); 422 } 423 PetscCall(MatMultTransposeAddEqual(A,S,10,&flg)); 424 if (!flg) { 425 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error axpy mult add (T)\n",test)); 426 } 427 PetscCall(MatConvert(A,MATDENSE,MAT_INITIAL_MATRIX,&dA)); 428 PetscCall(MatConvert(S,MATDENSE,MAT_INITIAL_MATRIX,&dS)); 429 PetscCall(MatAXPY(dA,-1.0,dS,SAME_NONZERO_PATTERN)); 430 PetscCall(MatNorm(dA,NORM_FROBENIUS,&err)); 431 if (err >= tol) { 432 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error mat submatrix %g\n",test,(double)err)); 433 } 434 PetscCall(MatDestroy(&dA)); 435 PetscCall(MatDestroy(&dS)); 436 } 437 438 if (testreset && (ntest == 1 || test == ntest-2)) { 439 /* reset MATSHELL */ 440 PetscCall(MatAssemblyBegin(S,MAT_FINAL_ASSEMBLY)); 441 PetscCall(MatAssemblyEnd(S,MAT_FINAL_ASSEMBLY)); 442 /* reset A */ 443 PetscCall(MatCopy(user->B,A,DIFFERENT_NONZERO_PATTERN)); 444 } 445 } 446 447 PetscCall(MatDestroy(&A)); 448 PetscCall(MatDestroy(&S)); 449 PetscCall(PetscFree(user)); 450 PetscCall(PetscFinalize()); 451 return 0; 452 } 453 454 /*TEST 455 456 testset: 457 suffix: rect 458 requires: !single 459 output_file: output/ex221_1.out 460 nsize: {{1 3}} 461 args: -loop 3 -keep {{0 1}} -M {{12 19}} -N {{19 12}} -submat {{0 1}} -test_axpy_different {{0 1}} 462 463 testset: 464 suffix: square 465 requires: !single 466 output_file: output/ex221_1.out 467 nsize: {{1 3}} 468 args: -M 21 -N 21 -loop 4 -rows_only {{0 1}} -keep {{0 1}} -submat {{0 1}} -test_axpy_different {{0 1}} 469 TEST*/ 470