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 CHKERRQ(MatShellGetContext(A,&user)); 16 CHKERRQ(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 CHKERRQ(MatShellGetContext(A,&user)); 26 CHKERRQ(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 CHKERRQ(MatShellGetContext(A,&user)); 36 CHKERRQ(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 CHKERRQ(MatShellGetContext(A,&user)); 46 CHKERRQ(MatShellGetContext(X,&userX)); 47 PetscCheckFalse(user != userX,PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"This should not happen"); 48 CHKERRQ(PetscObjectReference((PetscObject)user->B)); 49 PetscFunctionReturn(0); 50 } 51 52 static PetscErrorCode MatDestroy_User(Mat A) 53 { 54 User user; 55 56 PetscFunctionBegin; 57 CHKERRQ(MatShellGetContext(A,&user)); 58 CHKERRQ(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 CHKERRQ(PetscInitialize(&argc,&args,(char*)0,help)); 80 CHKERRMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 81 CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 82 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL)); 83 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL)); 84 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-ml",&m,NULL)); 85 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-nl",&n,NULL)); 86 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-square_nc",&nc,NULL)); 87 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-rows_only",&ronl,NULL)); 88 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-randomize",&randomize,NULL)); 89 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-submat",&submat,NULL)); 90 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-test_zerorows",&testzerorows,NULL)); 91 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-test_diagscale",&testdiagscale,NULL)); 92 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-test_getdiag",&testgetdiag,NULL)); 93 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-test_shift",&testshift,NULL)); 94 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-test_scale",&testscale,NULL)); 95 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-test_dup",&testdup,NULL)); 96 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-test_reset",&testreset,NULL)); 97 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-test_submat",&testsubmat,NULL)); 98 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-test_axpy",&testaxpy,NULL)); 99 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-test_axpy_different",&testaxpyd,NULL)); 100 CHKERRQ(PetscOptionsGetBool(NULL,NULL,"-test_axpy_error",&testaxpyerr,NULL)); 101 CHKERRQ(PetscOptionsGetInt(NULL,NULL,"-loop",&ntest,NULL)); 102 CHKERRQ(PetscOptionsGetReal(NULL,NULL,"-tol",&tol,NULL)); 103 CHKERRQ(PetscOptionsGetScalar(NULL,NULL,"-diag",&diag,NULL)); 104 CHKERRQ(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 CHKERRQ(MatCreateDense(PETSC_COMM_WORLD,m,n,M,N,NULL,&A)); 114 CHKERRQ(MatGetLocalSize(A,&m,&n)); 115 CHKERRQ(MatGetSize(A,&M,&N)); 116 CHKERRQ(MatGetOwnershipRange(A,&s1,NULL)); 117 s2 = 1; 118 while (s2 < M) s2 *= 10; 119 CHKERRQ(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 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 126 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 127 128 if (submat) { 129 Mat A2; 130 IS r,c; 131 PetscInt rst,ren,cst,cen; 132 133 CHKERRQ(MatGetOwnershipRange(A,&rst,&ren)); 134 CHKERRQ(MatGetOwnershipRangeColumn(A,&cst,&cen)); 135 CHKERRQ(ISCreateStride(PetscObjectComm((PetscObject)A),(ren-rst)/2,rst,1,&r)); 136 CHKERRQ(ISCreateStride(PetscObjectComm((PetscObject)A),(cen-cst)/2,cst,1,&c)); 137 CHKERRQ(MatCreateSubMatrix(A,r,c,MAT_INITIAL_MATRIX,&A2)); 138 CHKERRQ(ISDestroy(&r)); 139 CHKERRQ(ISDestroy(&c)); 140 CHKERRQ(MatDestroy(&A)); 141 A = A2; 142 } 143 144 CHKERRQ(MatGetSize(A,&M,&N)); 145 CHKERRQ(MatGetLocalSize(A,&m,&n)); 146 CHKERRQ(MatHasCongruentLayouts(A,&cong)); 147 148 CHKERRQ(MatConvert(A,MATAIJ,MAT_INPLACE_MATRIX,&A)); 149 CHKERRQ(MatSetOption(A,MAT_KEEP_NONZERO_PATTERN,keep)); 150 CHKERRQ(PetscObjectSetName((PetscObject)A,"initial")); 151 CHKERRQ(MatViewFromOptions(A,NULL,"-view_mat")); 152 153 CHKERRQ(PetscNew(&user)); 154 CHKERRQ(MatCreateShell(PETSC_COMM_WORLD,m,n,M,N,user,&S)); 155 CHKERRQ(MatShellSetOperation(S,MATOP_MULT,(void (*)(void))MatMult_User)); 156 CHKERRQ(MatShellSetOperation(S,MATOP_MULT_TRANSPOSE,(void (*)(void))MatMultTranspose_User)); 157 if (cong) { 158 CHKERRQ(MatShellSetOperation(S,MATOP_GET_DIAGONAL,(void (*)(void))MatGetDiagonal_User)); 159 } 160 CHKERRQ(MatShellSetOperation(S,MATOP_COPY,(void (*)(void))MatCopy_User)); 161 CHKERRQ(MatShellSetOperation(S,MATOP_DESTROY,(void (*)(void))MatDestroy_User)); 162 CHKERRQ(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 CHKERRQ(MatMultAddEqual(A,S,10,&flg)); 171 if (!flg) { 172 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error mult add\n",test)); 173 } 174 CHKERRQ(MatMultTransposeAddEqual(A,S,10,&flg)); 175 if (!flg) { 176 CHKERRQ(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 CHKERRQ(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 CHKERRQ(ISCreateGeneral(PETSC_COMM_WORLD,nr,idxs,PETSC_OWN_POINTER,&zr)); 195 CHKERRQ(PetscObjectSetName((PetscObject)zr,"ZR")); 196 CHKERRQ(ISViewFromOptions(zr,NULL,"-view_is")); 197 CHKERRQ(MatCreateVecs(A,&x,&b1)); 198 if (randomize) { 199 CHKERRQ(VecSetRandom(x,NULL)); 200 CHKERRQ(VecSetRandom(b1,NULL)); 201 } else { 202 CHKERRQ(VecSet(x,11.4)); 203 CHKERRQ(VecSet(b1,-14.2)); 204 } 205 CHKERRQ(VecDuplicate(b1,&b2)); 206 CHKERRQ(VecCopy(b1,b2)); 207 CHKERRQ(PetscObjectSetName((PetscObject)b1,"A_B1")); 208 CHKERRQ(PetscObjectSetName((PetscObject)b2,"A_B2")); 209 if (size > 1 && !cong) { /* MATMPIAIJ ZeroRows and ZeroRowsColumns are buggy in this case */ 210 CHKERRQ(VecDestroy(&b1)); 211 } 212 if (ronl) { 213 CHKERRQ(MatZeroRowsIS(A,zr,diag,x,b1)); 214 CHKERRQ(MatZeroRowsIS(S,zr,diag,x,b2)); 215 } else { 216 CHKERRQ(MatZeroRowsColumnsIS(A,zr,diag,x,b1)); 217 CHKERRQ(MatZeroRowsColumnsIS(S,zr,diag,x,b2)); 218 CHKERRQ(ISDestroy(&zr)); 219 /* Mix zerorows and zerorowscols */ 220 nr = 0; 221 idxs = NULL; 222 if (rank == 0) { 223 nr = 1; 224 CHKERRQ(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 CHKERRQ(ISCreateGeneral(PETSC_COMM_WORLD,nr,idxs,PETSC_OWN_POINTER,&zr)); 233 CHKERRQ(PetscObjectSetName((PetscObject)zr,"ZR2")); 234 CHKERRQ(ISViewFromOptions(zr,NULL,"-view_is")); 235 CHKERRQ(MatZeroRowsIS(A,zr,diag*2.0+PETSC_SMALL,NULL,NULL)); 236 CHKERRQ(MatZeroRowsIS(S,zr,diag*2.0+PETSC_SMALL,NULL,NULL)); 237 } 238 CHKERRQ(ISDestroy(&zr)); 239 240 if (b1) { 241 Vec b; 242 243 CHKERRQ(VecViewFromOptions(b1,NULL,"-view_b")); 244 CHKERRQ(VecViewFromOptions(b2,NULL,"-view_b")); 245 CHKERRQ(VecDuplicate(b1,&b)); 246 CHKERRQ(VecCopy(b1,b)); 247 CHKERRQ(VecAXPY(b,-1.0,b2)); 248 CHKERRQ(VecNorm(b,NORM_INFINITY,&err)); 249 if (err >= tol) { 250 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error b %g\n",test,(double)err)); 251 } 252 CHKERRQ(VecDestroy(&b)); 253 } 254 CHKERRQ(VecDestroy(&b1)); 255 CHKERRQ(VecDestroy(&b2)); 256 CHKERRQ(VecDestroy(&x)); 257 CHKERRQ(MatConvert(S,MATDENSE,MAT_INITIAL_MATRIX,&B)); 258 259 CHKERRQ(MatCreateTranspose(S,&ST)); 260 CHKERRQ(MatComputeOperator(ST,MATDENSE,&BT)); 261 CHKERRQ(MatTranspose(BT,MAT_INITIAL_MATRIX,&BTT)); 262 CHKERRQ(PetscObjectSetName((PetscObject)B,"S")); 263 CHKERRQ(PetscObjectSetName((PetscObject)BTT,"STT")); 264 CHKERRQ(MatConvert(A,MATDENSE,MAT_INITIAL_MATRIX,&C)); 265 CHKERRQ(PetscObjectSetName((PetscObject)C,"A")); 266 267 CHKERRQ(MatViewFromOptions(C,NULL,"-view_mat")); 268 CHKERRQ(MatViewFromOptions(B,NULL,"-view_mat")); 269 CHKERRQ(MatViewFromOptions(BTT,NULL,"-view_mat")); 270 271 CHKERRQ(MatAXPY(C,-1.0,B,SAME_NONZERO_PATTERN)); 272 CHKERRQ(MatNorm(C,NORM_FROBENIUS,&err)); 273 if (err >= tol) { 274 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error mat mult after %s %g\n",test,ronl ? "MatZeroRows" : "MatZeroRowsColumns",(double)err)); 275 } 276 277 CHKERRQ(MatConvert(A,MATDENSE,MAT_REUSE_MATRIX,&C)); 278 CHKERRQ(MatAXPY(C,-1.0,BTT,SAME_NONZERO_PATTERN)); 279 CHKERRQ(MatNorm(C,NORM_FROBENIUS,&err)); 280 if (err >= tol) { 281 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error mat mult transpose after %s %g\n",test,ronl ? "MatZeroRows" : "MatZeroRowsColumns",(double)err)); 282 } 283 284 CHKERRQ(MatDestroy(&ST)); 285 CHKERRQ(MatDestroy(&BTT)); 286 CHKERRQ(MatDestroy(&BT)); 287 CHKERRQ(MatDestroy(&B)); 288 CHKERRQ(MatDestroy(&C)); 289 } 290 if (testdiagscale) { /* MatDiagonalScale() */ 291 Vec vr,vl; 292 293 CHKERRQ(MatCreateVecs(A,&vr,&vl)); 294 if (randomize) { 295 CHKERRQ(VecSetRandom(vr,NULL)); 296 CHKERRQ(VecSetRandom(vl,NULL)); 297 } else { 298 CHKERRQ(VecSet(vr,test%2 ? 0.15 : 1.0/0.15)); 299 CHKERRQ(VecSet(vl,test%2 ? -1.2 : 1.0/-1.2)); 300 } 301 CHKERRQ(MatDiagonalScale(A,vl,vr)); 302 CHKERRQ(MatDiagonalScale(S,vl,vr)); 303 CHKERRQ(VecDestroy(&vr)); 304 CHKERRQ(VecDestroy(&vl)); 305 } 306 307 if (testscale) { /* MatScale() */ 308 CHKERRQ(MatScale(A,test%2 ? 1.4 : 1.0/1.4)); 309 CHKERRQ(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 CHKERRQ(MatShift(A,test%2 ? -77.5 : 77.5)); 314 CHKERRQ(MatShift(S,test%2 ? -77.5 : 77.5)); 315 } 316 317 if (testgetdiag && cong) { /* MatGetDiagonal() */ 318 Vec dA,dS; 319 320 CHKERRQ(MatCreateVecs(A,&dA,NULL)); 321 CHKERRQ(MatCreateVecs(S,&dS,NULL)); 322 CHKERRQ(MatGetDiagonal(A,dA)); 323 CHKERRQ(MatGetDiagonal(S,dS)); 324 CHKERRQ(VecAXPY(dA,-1.0,dS)); 325 CHKERRQ(VecNorm(dA,NORM_INFINITY,&err)); 326 if (err >= tol) { 327 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error diag %g\n",test,(double)err)); 328 } 329 CHKERRQ(VecDestroy(&dA)); 330 CHKERRQ(VecDestroy(&dS)); 331 } 332 333 if (testdup && !test) { 334 Mat A2, S2; 335 336 CHKERRQ(MatDuplicate(A,MAT_COPY_VALUES,&A2)); 337 CHKERRQ(MatDuplicate(S,MAT_COPY_VALUES,&S2)); 338 CHKERRQ(MatDestroy(&A)); 339 CHKERRQ(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 CHKERRQ(MatGetOwnershipRange(A,&rst,&ren)); 350 CHKERRQ(MatGetOwnershipRangeColumn(A,&cst,&cen)); 351 CHKERRQ(ISCreateStride(PetscObjectComm((PetscObject)A),(ren-rst)/2,rst,1,&r)); 352 CHKERRQ(ISCreateStride(PetscObjectComm((PetscObject)A),(cen-cst)/2,cst,1,&c)); 353 CHKERRQ(MatCreateSubMatrix(A,r,c,MAT_INITIAL_MATRIX,&sA)); 354 CHKERRQ(MatCreateSubMatrix(S,r,c,MAT_INITIAL_MATRIX,&sS)); 355 CHKERRQ(MatMultAddEqual(sA,sS,10,&flg)); 356 if (!flg) { 357 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error submatrix mult add\n",test)); 358 } 359 CHKERRQ(MatMultTransposeAddEqual(sA,sS,10,&flg)); 360 if (!flg) { 361 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error submatrix mult add (T)\n",test)); 362 } 363 CHKERRQ(MatConvert(sA,MATDENSE,MAT_INITIAL_MATRIX,&dA)); 364 CHKERRQ(MatConvert(sS,MATDENSE,MAT_INITIAL_MATRIX,&dS)); 365 CHKERRQ(MatAXPY(dA,-1.0,dS,SAME_NONZERO_PATTERN)); 366 CHKERRQ(MatNorm(dA,NORM_FROBENIUS,&err)); 367 if (err >= tol) { 368 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error mat submatrix %g\n",test,(double)err)); 369 } 370 CHKERRQ(MatDestroy(&sA)); 371 CHKERRQ(MatDestroy(&sS)); 372 CHKERRQ(MatDestroy(&dA)); 373 CHKERRQ(MatDestroy(&dS)); 374 CHKERRQ(MatCreateTranspose(A,&At)); 375 CHKERRQ(MatCreateTranspose(S,&St)); 376 CHKERRQ(MatCreateSubMatrix(At,c,r,MAT_INITIAL_MATRIX,&sA)); 377 CHKERRQ(MatCreateSubMatrix(St,c,r,MAT_INITIAL_MATRIX,&sS)); 378 CHKERRQ(MatMultAddEqual(sA,sS,10,&flg)); 379 if (!flg) { 380 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error submatrix (T) mult add\n",test)); 381 } 382 CHKERRQ(MatMultTransposeAddEqual(sA,sS,10,&flg)); 383 if (!flg) { 384 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error submatrix (T) mult add (T)\n",test)); 385 } 386 CHKERRQ(MatConvert(sA,MATDENSE,MAT_INITIAL_MATRIX,&dA)); 387 CHKERRQ(MatConvert(sS,MATDENSE,MAT_INITIAL_MATRIX,&dS)); 388 CHKERRQ(MatAXPY(dA,-1.0,dS,SAME_NONZERO_PATTERN)); 389 CHKERRQ(MatNorm(dA,NORM_FROBENIUS,&err)); 390 if (err >= tol) { 391 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error mat submatrix (T) %g\n",test,(double)err)); 392 } 393 CHKERRQ(MatDestroy(&sA)); 394 CHKERRQ(MatDestroy(&sS)); 395 CHKERRQ(MatDestroy(&dA)); 396 CHKERRQ(MatDestroy(&dS)); 397 CHKERRQ(MatDestroy(&At)); 398 CHKERRQ(MatDestroy(&St)); 399 CHKERRQ(ISDestroy(&r)); 400 CHKERRQ(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 CHKERRQ(MatDuplicate(A,MAT_COPY_VALUES,&tA)); 408 if (testaxpyd && !(test%2)) { 409 CHKERRQ(PetscObjectReference((PetscObject)tA)); 410 tS = tA; 411 } else { 412 CHKERRQ(PetscObjectReference((PetscObject)S)); 413 tS = S; 414 } 415 CHKERRQ(MatAXPY(A,0.5,tA,str[test%3])); 416 CHKERRQ(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) CHKERRQ(MatScale(tA,0)); 419 CHKERRQ(MatDestroy(&tA)); 420 CHKERRQ(MatDestroy(&tS)); 421 CHKERRQ(MatMultAddEqual(A,S,10,&flg)); 422 if (!flg) { 423 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error axpy mult add\n",test)); 424 } 425 CHKERRQ(MatMultTransposeAddEqual(A,S,10,&flg)); 426 if (!flg) { 427 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error axpy mult add (T)\n",test)); 428 } 429 CHKERRQ(MatConvert(A,MATDENSE,MAT_INITIAL_MATRIX,&dA)); 430 CHKERRQ(MatConvert(S,MATDENSE,MAT_INITIAL_MATRIX,&dS)); 431 CHKERRQ(MatAXPY(dA,-1.0,dS,SAME_NONZERO_PATTERN)); 432 CHKERRQ(MatNorm(dA,NORM_FROBENIUS,&err)); 433 if (err >= tol) { 434 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error mat submatrix %g\n",test,(double)err)); 435 } 436 CHKERRQ(MatDestroy(&dA)); 437 CHKERRQ(MatDestroy(&dS)); 438 } 439 440 if (testreset && (ntest == 1 || test == ntest-2)) { 441 /* reset MATSHELL */ 442 CHKERRQ(MatAssemblyBegin(S,MAT_FINAL_ASSEMBLY)); 443 CHKERRQ(MatAssemblyEnd(S,MAT_FINAL_ASSEMBLY)); 444 /* reset A */ 445 CHKERRQ(MatCopy(user->B,A,DIFFERENT_NONZERO_PATTERN)); 446 } 447 } 448 449 CHKERRQ(MatDestroy(&A)); 450 CHKERRQ(MatDestroy(&S)); 451 CHKERRQ(PetscFree(user)); 452 CHKERRQ(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