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