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