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