1c4762a1bSJed Brown static char help[] = "Tests various routines for MATSHELL\n\n"; 2c4762a1bSJed Brown 3c4762a1bSJed Brown #include <petscmat.h> 4c4762a1bSJed Brown 5c4762a1bSJed Brown typedef struct _n_User *User; 6c4762a1bSJed Brown struct _n_User { 7c4762a1bSJed Brown Mat B; 8c4762a1bSJed Brown }; 9c4762a1bSJed Brown 10c4762a1bSJed Brown static PetscErrorCode MatGetDiagonal_User(Mat A,Vec X) 11c4762a1bSJed Brown { 12c4762a1bSJed Brown User user; 13c4762a1bSJed Brown 14c4762a1bSJed Brown PetscFunctionBegin; 159566063dSJacob Faibussowitsch PetscCall(MatShellGetContext(A,&user)); 169566063dSJacob Faibussowitsch PetscCall(MatGetDiagonal(user->B,X)); 17c4762a1bSJed Brown PetscFunctionReturn(0); 18c4762a1bSJed Brown } 19c4762a1bSJed Brown 20c4762a1bSJed Brown static PetscErrorCode MatMult_User(Mat A,Vec X,Vec Y) 21c4762a1bSJed Brown { 22c4762a1bSJed Brown User user; 23c4762a1bSJed Brown 24c4762a1bSJed Brown PetscFunctionBegin; 259566063dSJacob Faibussowitsch PetscCall(MatShellGetContext(A,&user)); 269566063dSJacob Faibussowitsch PetscCall(MatMult(user->B,X,Y)); 27c4762a1bSJed Brown PetscFunctionReturn(0); 28c4762a1bSJed Brown } 29c4762a1bSJed Brown 30c4762a1bSJed Brown static PetscErrorCode MatMultTranspose_User(Mat A,Vec X,Vec Y) 31c4762a1bSJed Brown { 32c4762a1bSJed Brown User user; 33c4762a1bSJed Brown 34c4762a1bSJed Brown PetscFunctionBegin; 359566063dSJacob Faibussowitsch PetscCall(MatShellGetContext(A,&user)); 369566063dSJacob Faibussowitsch PetscCall(MatMultTranspose(user->B,X,Y)); 37c4762a1bSJed Brown PetscFunctionReturn(0); 38c4762a1bSJed Brown } 39c4762a1bSJed Brown 40c4762a1bSJed Brown static PetscErrorCode MatCopy_User(Mat A,Mat X,MatStructure str) 41c4762a1bSJed Brown { 42c4762a1bSJed Brown User user,userX; 43c4762a1bSJed Brown 44c4762a1bSJed Brown PetscFunctionBegin; 459566063dSJacob Faibussowitsch PetscCall(MatShellGetContext(A,&user)); 469566063dSJacob Faibussowitsch PetscCall(MatShellGetContext(X,&userX)); 47*08401ef6SPierre Jolivet PetscCheck(user == userX,PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"This should not happen"); 489566063dSJacob Faibussowitsch PetscCall(PetscObjectReference((PetscObject)user->B)); 49c4762a1bSJed Brown PetscFunctionReturn(0); 50c4762a1bSJed Brown } 51c4762a1bSJed Brown 52c4762a1bSJed Brown static PetscErrorCode MatDestroy_User(Mat A) 53c4762a1bSJed Brown { 54c4762a1bSJed Brown User user; 55c4762a1bSJed Brown 56c4762a1bSJed Brown PetscFunctionBegin; 579566063dSJacob Faibussowitsch PetscCall(MatShellGetContext(A,&user)); 589566063dSJacob Faibussowitsch PetscCall(PetscObjectDereference((PetscObject)user->B)); 59c4762a1bSJed Brown PetscFunctionReturn(0); 60c4762a1bSJed Brown } 61c4762a1bSJed Brown 62c4762a1bSJed Brown int main(int argc,char **args) 63c4762a1bSJed Brown { 64c4762a1bSJed Brown User user; 65c4762a1bSJed Brown Mat A,S; 66c4762a1bSJed Brown PetscScalar *data,diag = 1.3; 67c4762a1bSJed Brown PetscReal tol = PETSC_SMALL; 68c4762a1bSJed Brown PetscInt i,j,m = PETSC_DECIDE,n = PETSC_DECIDE,M = 17,N = 15,s1,s2; 69c4762a1bSJed Brown PetscInt test, ntest = 2; 70c4762a1bSJed Brown PetscMPIInt rank,size; 71c4762a1bSJed Brown PetscBool nc = PETSC_FALSE, cong, flg; 72c4762a1bSJed Brown PetscBool ronl = PETSC_TRUE; 73c4762a1bSJed Brown PetscBool randomize = PETSC_FALSE, submat = PETSC_FALSE; 74c4762a1bSJed Brown PetscBool keep = PETSC_FALSE; 75c4762a1bSJed Brown PetscBool testzerorows = PETSC_TRUE, testdiagscale = PETSC_TRUE, testgetdiag = PETSC_TRUE, testsubmat = PETSC_TRUE; 76c4762a1bSJed Brown PetscBool testshift = PETSC_TRUE, testscale = PETSC_TRUE, testdup = PETSC_TRUE, testreset = PETSC_TRUE; 77c4762a1bSJed Brown PetscBool testaxpy = PETSC_TRUE, testaxpyd = PETSC_TRUE, testaxpyerr = PETSC_FALSE; 78c4762a1bSJed Brown 799566063dSJacob Faibussowitsch PetscCall(PetscInitialize(&argc,&args,(char*)0,help)); 809566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 819566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 829566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL)); 839566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL)); 849566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetInt(NULL,NULL,"-ml",&m,NULL)); 859566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetInt(NULL,NULL,"-nl",&n,NULL)); 869566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL,NULL,"-square_nc",&nc,NULL)); 879566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL,NULL,"-rows_only",&ronl,NULL)); 889566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL,NULL,"-randomize",&randomize,NULL)); 899566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL,NULL,"-submat",&submat,NULL)); 909566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_zerorows",&testzerorows,NULL)); 919566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_diagscale",&testdiagscale,NULL)); 929566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_getdiag",&testgetdiag,NULL)); 939566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_shift",&testshift,NULL)); 949566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_scale",&testscale,NULL)); 959566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_dup",&testdup,NULL)); 969566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_reset",&testreset,NULL)); 979566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_submat",&testsubmat,NULL)); 989566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_axpy",&testaxpy,NULL)); 999566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_axpy_different",&testaxpyd,NULL)); 1009566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_axpy_error",&testaxpyerr,NULL)); 1019566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetInt(NULL,NULL,"-loop",&ntest,NULL)); 1029566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetReal(NULL,NULL,"-tol",&tol,NULL)); 1039566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetScalar(NULL,NULL,"-diag",&diag,NULL)); 1049566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL,NULL,"-keep",&keep,NULL)); 105c4762a1bSJed Brown /* This tests square matrices with different row/col layout */ 106c4762a1bSJed Brown if (nc && size > 1) { 107c4762a1bSJed Brown M = PetscMax(PetscMax(N,M),1); 108c4762a1bSJed Brown N = M; 109c4762a1bSJed Brown m = n = 0; 110c4762a1bSJed Brown if (rank == 0) { m = M-1; n = 1; } 111c4762a1bSJed Brown else if (rank == 1) { m = 1; n = N-1; } 112c4762a1bSJed Brown } 1139566063dSJacob Faibussowitsch PetscCall(MatCreateDense(PETSC_COMM_WORLD,m,n,M,N,NULL,&A)); 1149566063dSJacob Faibussowitsch PetscCall(MatGetLocalSize(A,&m,&n)); 1159566063dSJacob Faibussowitsch PetscCall(MatGetSize(A,&M,&N)); 1169566063dSJacob Faibussowitsch PetscCall(MatGetOwnershipRange(A,&s1,NULL)); 117c4762a1bSJed Brown s2 = 1; 118c4762a1bSJed Brown while (s2 < M) s2 *= 10; 1199566063dSJacob Faibussowitsch PetscCall(MatDenseGetArray(A,&data)); 120c4762a1bSJed Brown for (j = 0; j < N; j++) { 121c4762a1bSJed Brown for (i = 0; i < m; i++) { 122c4762a1bSJed Brown data[j*m + i] = s2*j + i + s1 + 1; 123c4762a1bSJed Brown } 124c4762a1bSJed Brown } 1259566063dSJacob Faibussowitsch PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 1269566063dSJacob Faibussowitsch PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 127c4762a1bSJed Brown 128c4762a1bSJed Brown if (submat) { 129c4762a1bSJed Brown Mat A2; 130c4762a1bSJed Brown IS r,c; 131c4762a1bSJed Brown PetscInt rst,ren,cst,cen; 132c4762a1bSJed Brown 1339566063dSJacob Faibussowitsch PetscCall(MatGetOwnershipRange(A,&rst,&ren)); 1349566063dSJacob Faibussowitsch PetscCall(MatGetOwnershipRangeColumn(A,&cst,&cen)); 1359566063dSJacob Faibussowitsch PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A),(ren-rst)/2,rst,1,&r)); 1369566063dSJacob Faibussowitsch PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A),(cen-cst)/2,cst,1,&c)); 1379566063dSJacob Faibussowitsch PetscCall(MatCreateSubMatrix(A,r,c,MAT_INITIAL_MATRIX,&A2)); 1389566063dSJacob Faibussowitsch PetscCall(ISDestroy(&r)); 1399566063dSJacob Faibussowitsch PetscCall(ISDestroy(&c)); 1409566063dSJacob Faibussowitsch PetscCall(MatDestroy(&A)); 141c4762a1bSJed Brown A = A2; 142c4762a1bSJed Brown } 143c4762a1bSJed Brown 1449566063dSJacob Faibussowitsch PetscCall(MatGetSize(A,&M,&N)); 1459566063dSJacob Faibussowitsch PetscCall(MatGetLocalSize(A,&m,&n)); 1469566063dSJacob Faibussowitsch PetscCall(MatHasCongruentLayouts(A,&cong)); 147c4762a1bSJed Brown 1489566063dSJacob Faibussowitsch PetscCall(MatConvert(A,MATAIJ,MAT_INPLACE_MATRIX,&A)); 1499566063dSJacob Faibussowitsch PetscCall(MatSetOption(A,MAT_KEEP_NONZERO_PATTERN,keep)); 1509566063dSJacob Faibussowitsch PetscCall(PetscObjectSetName((PetscObject)A,"initial")); 1519566063dSJacob Faibussowitsch PetscCall(MatViewFromOptions(A,NULL,"-view_mat")); 152c4762a1bSJed Brown 1539566063dSJacob Faibussowitsch PetscCall(PetscNew(&user)); 1549566063dSJacob Faibussowitsch PetscCall(MatCreateShell(PETSC_COMM_WORLD,m,n,M,N,user,&S)); 1559566063dSJacob Faibussowitsch PetscCall(MatShellSetOperation(S,MATOP_MULT,(void (*)(void))MatMult_User)); 1569566063dSJacob Faibussowitsch PetscCall(MatShellSetOperation(S,MATOP_MULT_TRANSPOSE,(void (*)(void))MatMultTranspose_User)); 157c4762a1bSJed Brown if (cong) { 1589566063dSJacob Faibussowitsch PetscCall(MatShellSetOperation(S,MATOP_GET_DIAGONAL,(void (*)(void))MatGetDiagonal_User)); 159c4762a1bSJed Brown } 1609566063dSJacob Faibussowitsch PetscCall(MatShellSetOperation(S,MATOP_COPY,(void (*)(void))MatCopy_User)); 1619566063dSJacob Faibussowitsch PetscCall(MatShellSetOperation(S,MATOP_DESTROY,(void (*)(void))MatDestroy_User)); 1629566063dSJacob Faibussowitsch PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&user->B)); 163c4762a1bSJed Brown 164c4762a1bSJed Brown /* Square and rows only scaling */ 165c4762a1bSJed Brown ronl = cong ? ronl : PETSC_TRUE; 166c4762a1bSJed Brown 167c4762a1bSJed Brown for (test = 0; test < ntest; test++) { 168c4762a1bSJed Brown PetscReal err; 169c4762a1bSJed Brown 1709566063dSJacob Faibussowitsch PetscCall(MatMultAddEqual(A,S,10,&flg)); 171c4762a1bSJed Brown if (!flg) { 1729566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error mult add\n",test)); 173c4762a1bSJed Brown } 1749566063dSJacob Faibussowitsch PetscCall(MatMultTransposeAddEqual(A,S,10,&flg)); 175c4762a1bSJed Brown if (!flg) { 1769566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error mult add (T)\n",test)); 177c4762a1bSJed Brown } 178c4762a1bSJed Brown if (testzerorows) { 179c4762a1bSJed Brown Mat ST,B,C,BT,BTT; 180c4762a1bSJed Brown IS zr; 181c4762a1bSJed Brown Vec x = NULL, b1 = NULL, b2 = NULL; 182c4762a1bSJed Brown PetscInt *idxs = NULL, nr = 0; 183c4762a1bSJed Brown 184c4762a1bSJed Brown if (rank == (test%size)) { 185c4762a1bSJed Brown nr = 1; 1869566063dSJacob Faibussowitsch PetscCall(PetscMalloc1(nr,&idxs)); 187c4762a1bSJed Brown if (test%2) { 188c4762a1bSJed Brown idxs[0] = (2*M - 1 - test/2)%M; 189c4762a1bSJed Brown } else { 190c4762a1bSJed Brown idxs[0] = (test/2)%M; 191c4762a1bSJed Brown } 192c4762a1bSJed Brown idxs[0] = PetscMax(idxs[0],0); 193c4762a1bSJed Brown } 1949566063dSJacob Faibussowitsch PetscCall(ISCreateGeneral(PETSC_COMM_WORLD,nr,idxs,PETSC_OWN_POINTER,&zr)); 1959566063dSJacob Faibussowitsch PetscCall(PetscObjectSetName((PetscObject)zr,"ZR")); 1969566063dSJacob Faibussowitsch PetscCall(ISViewFromOptions(zr,NULL,"-view_is")); 1979566063dSJacob Faibussowitsch PetscCall(MatCreateVecs(A,&x,&b1)); 198c4762a1bSJed Brown if (randomize) { 1999566063dSJacob Faibussowitsch PetscCall(VecSetRandom(x,NULL)); 2009566063dSJacob Faibussowitsch PetscCall(VecSetRandom(b1,NULL)); 201c4762a1bSJed Brown } else { 2029566063dSJacob Faibussowitsch PetscCall(VecSet(x,11.4)); 2039566063dSJacob Faibussowitsch PetscCall(VecSet(b1,-14.2)); 204c4762a1bSJed Brown } 2059566063dSJacob Faibussowitsch PetscCall(VecDuplicate(b1,&b2)); 2069566063dSJacob Faibussowitsch PetscCall(VecCopy(b1,b2)); 2079566063dSJacob Faibussowitsch PetscCall(PetscObjectSetName((PetscObject)b1,"A_B1")); 2089566063dSJacob Faibussowitsch PetscCall(PetscObjectSetName((PetscObject)b2,"A_B2")); 209c4762a1bSJed Brown if (size > 1 && !cong) { /* MATMPIAIJ ZeroRows and ZeroRowsColumns are buggy in this case */ 2109566063dSJacob Faibussowitsch PetscCall(VecDestroy(&b1)); 211c4762a1bSJed Brown } 212c4762a1bSJed Brown if (ronl) { 2139566063dSJacob Faibussowitsch PetscCall(MatZeroRowsIS(A,zr,diag,x,b1)); 2149566063dSJacob Faibussowitsch PetscCall(MatZeroRowsIS(S,zr,diag,x,b2)); 215c4762a1bSJed Brown } else { 2169566063dSJacob Faibussowitsch PetscCall(MatZeroRowsColumnsIS(A,zr,diag,x,b1)); 2179566063dSJacob Faibussowitsch PetscCall(MatZeroRowsColumnsIS(S,zr,diag,x,b2)); 2189566063dSJacob Faibussowitsch PetscCall(ISDestroy(&zr)); 219c4762a1bSJed Brown /* Mix zerorows and zerorowscols */ 220c4762a1bSJed Brown nr = 0; 221c4762a1bSJed Brown idxs = NULL; 222dd400576SPatrick Sanan if (rank == 0) { 223c4762a1bSJed Brown nr = 1; 2249566063dSJacob Faibussowitsch PetscCall(PetscMalloc1(nr,&idxs)); 225c4762a1bSJed Brown if (test%2) { 226c4762a1bSJed Brown idxs[0] = (3*M - 2 - test/2)%M; 227c4762a1bSJed Brown } else { 228c4762a1bSJed Brown idxs[0] = (test/2+1)%M; 229c4762a1bSJed Brown } 230c4762a1bSJed Brown idxs[0] = PetscMax(idxs[0],0); 231c4762a1bSJed Brown } 2329566063dSJacob Faibussowitsch PetscCall(ISCreateGeneral(PETSC_COMM_WORLD,nr,idxs,PETSC_OWN_POINTER,&zr)); 2339566063dSJacob Faibussowitsch PetscCall(PetscObjectSetName((PetscObject)zr,"ZR2")); 2349566063dSJacob Faibussowitsch PetscCall(ISViewFromOptions(zr,NULL,"-view_is")); 2359566063dSJacob Faibussowitsch PetscCall(MatZeroRowsIS(A,zr,diag*2.0+PETSC_SMALL,NULL,NULL)); 2369566063dSJacob Faibussowitsch PetscCall(MatZeroRowsIS(S,zr,diag*2.0+PETSC_SMALL,NULL,NULL)); 237c4762a1bSJed Brown } 2389566063dSJacob Faibussowitsch PetscCall(ISDestroy(&zr)); 239c4762a1bSJed Brown 240c4762a1bSJed Brown if (b1) { 241c4762a1bSJed Brown Vec b; 242c4762a1bSJed Brown 2439566063dSJacob Faibussowitsch PetscCall(VecViewFromOptions(b1,NULL,"-view_b")); 2449566063dSJacob Faibussowitsch PetscCall(VecViewFromOptions(b2,NULL,"-view_b")); 2459566063dSJacob Faibussowitsch PetscCall(VecDuplicate(b1,&b)); 2469566063dSJacob Faibussowitsch PetscCall(VecCopy(b1,b)); 2479566063dSJacob Faibussowitsch PetscCall(VecAXPY(b,-1.0,b2)); 2489566063dSJacob Faibussowitsch PetscCall(VecNorm(b,NORM_INFINITY,&err)); 249c4762a1bSJed Brown if (err >= tol) { 2509566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error b %g\n",test,(double)err)); 251c4762a1bSJed Brown } 2529566063dSJacob Faibussowitsch PetscCall(VecDestroy(&b)); 253c4762a1bSJed Brown } 2549566063dSJacob Faibussowitsch PetscCall(VecDestroy(&b1)); 2559566063dSJacob Faibussowitsch PetscCall(VecDestroy(&b2)); 2569566063dSJacob Faibussowitsch PetscCall(VecDestroy(&x)); 2579566063dSJacob Faibussowitsch PetscCall(MatConvert(S,MATDENSE,MAT_INITIAL_MATRIX,&B)); 258c4762a1bSJed Brown 2599566063dSJacob Faibussowitsch PetscCall(MatCreateTranspose(S,&ST)); 2609566063dSJacob Faibussowitsch PetscCall(MatComputeOperator(ST,MATDENSE,&BT)); 2619566063dSJacob Faibussowitsch PetscCall(MatTranspose(BT,MAT_INITIAL_MATRIX,&BTT)); 2629566063dSJacob Faibussowitsch PetscCall(PetscObjectSetName((PetscObject)B,"S")); 2639566063dSJacob Faibussowitsch PetscCall(PetscObjectSetName((PetscObject)BTT,"STT")); 2649566063dSJacob Faibussowitsch PetscCall(MatConvert(A,MATDENSE,MAT_INITIAL_MATRIX,&C)); 2659566063dSJacob Faibussowitsch PetscCall(PetscObjectSetName((PetscObject)C,"A")); 266c4762a1bSJed Brown 2679566063dSJacob Faibussowitsch PetscCall(MatViewFromOptions(C,NULL,"-view_mat")); 2689566063dSJacob Faibussowitsch PetscCall(MatViewFromOptions(B,NULL,"-view_mat")); 2699566063dSJacob Faibussowitsch PetscCall(MatViewFromOptions(BTT,NULL,"-view_mat")); 270c4762a1bSJed Brown 2719566063dSJacob Faibussowitsch PetscCall(MatAXPY(C,-1.0,B,SAME_NONZERO_PATTERN)); 2729566063dSJacob Faibussowitsch PetscCall(MatNorm(C,NORM_FROBENIUS,&err)); 273c4762a1bSJed Brown if (err >= tol) { 2749566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error mat mult after %s %g\n",test,ronl ? "MatZeroRows" : "MatZeroRowsColumns",(double)err)); 275c4762a1bSJed Brown } 276c4762a1bSJed Brown 2779566063dSJacob Faibussowitsch PetscCall(MatConvert(A,MATDENSE,MAT_REUSE_MATRIX,&C)); 2789566063dSJacob Faibussowitsch PetscCall(MatAXPY(C,-1.0,BTT,SAME_NONZERO_PATTERN)); 2799566063dSJacob Faibussowitsch PetscCall(MatNorm(C,NORM_FROBENIUS,&err)); 280c4762a1bSJed Brown if (err >= tol) { 2819566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error mat mult transpose after %s %g\n",test,ronl ? "MatZeroRows" : "MatZeroRowsColumns",(double)err)); 282c4762a1bSJed Brown } 283c4762a1bSJed Brown 2849566063dSJacob Faibussowitsch PetscCall(MatDestroy(&ST)); 2859566063dSJacob Faibussowitsch PetscCall(MatDestroy(&BTT)); 2869566063dSJacob Faibussowitsch PetscCall(MatDestroy(&BT)); 2879566063dSJacob Faibussowitsch PetscCall(MatDestroy(&B)); 2889566063dSJacob Faibussowitsch PetscCall(MatDestroy(&C)); 289c4762a1bSJed Brown } 290c4762a1bSJed Brown if (testdiagscale) { /* MatDiagonalScale() */ 291c4762a1bSJed Brown Vec vr,vl; 292c4762a1bSJed Brown 2939566063dSJacob Faibussowitsch PetscCall(MatCreateVecs(A,&vr,&vl)); 294c4762a1bSJed Brown if (randomize) { 2959566063dSJacob Faibussowitsch PetscCall(VecSetRandom(vr,NULL)); 2969566063dSJacob Faibussowitsch PetscCall(VecSetRandom(vl,NULL)); 297c4762a1bSJed Brown } else { 2989566063dSJacob Faibussowitsch PetscCall(VecSet(vr,test%2 ? 0.15 : 1.0/0.15)); 2999566063dSJacob Faibussowitsch PetscCall(VecSet(vl,test%2 ? -1.2 : 1.0/-1.2)); 300c4762a1bSJed Brown } 3019566063dSJacob Faibussowitsch PetscCall(MatDiagonalScale(A,vl,vr)); 3029566063dSJacob Faibussowitsch PetscCall(MatDiagonalScale(S,vl,vr)); 3039566063dSJacob Faibussowitsch PetscCall(VecDestroy(&vr)); 3049566063dSJacob Faibussowitsch PetscCall(VecDestroy(&vl)); 305c4762a1bSJed Brown } 306c4762a1bSJed Brown 307c4762a1bSJed Brown if (testscale) { /* MatScale() */ 3089566063dSJacob Faibussowitsch PetscCall(MatScale(A,test%2 ? 1.4 : 1.0/1.4)); 3099566063dSJacob Faibussowitsch PetscCall(MatScale(S,test%2 ? 1.4 : 1.0/1.4)); 310c4762a1bSJed Brown } 311c4762a1bSJed Brown 312c4762a1bSJed Brown if (testshift && cong) { /* MatShift() : MATSHELL shift is broken when row/cols layout are not congruent and left/right scaling have been applied */ 3139566063dSJacob Faibussowitsch PetscCall(MatShift(A,test%2 ? -77.5 : 77.5)); 3149566063dSJacob Faibussowitsch PetscCall(MatShift(S,test%2 ? -77.5 : 77.5)); 315c4762a1bSJed Brown } 316c4762a1bSJed Brown 317c4762a1bSJed Brown if (testgetdiag && cong) { /* MatGetDiagonal() */ 318c4762a1bSJed Brown Vec dA,dS; 319c4762a1bSJed Brown 3209566063dSJacob Faibussowitsch PetscCall(MatCreateVecs(A,&dA,NULL)); 3219566063dSJacob Faibussowitsch PetscCall(MatCreateVecs(S,&dS,NULL)); 3229566063dSJacob Faibussowitsch PetscCall(MatGetDiagonal(A,dA)); 3239566063dSJacob Faibussowitsch PetscCall(MatGetDiagonal(S,dS)); 3249566063dSJacob Faibussowitsch PetscCall(VecAXPY(dA,-1.0,dS)); 3259566063dSJacob Faibussowitsch PetscCall(VecNorm(dA,NORM_INFINITY,&err)); 326c4762a1bSJed Brown if (err >= tol) { 3279566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error diag %g\n",test,(double)err)); 328c4762a1bSJed Brown } 3299566063dSJacob Faibussowitsch PetscCall(VecDestroy(&dA)); 3309566063dSJacob Faibussowitsch PetscCall(VecDestroy(&dS)); 331c4762a1bSJed Brown } 332c4762a1bSJed Brown 333c4762a1bSJed Brown if (testdup && !test) { 334c4762a1bSJed Brown Mat A2, S2; 335c4762a1bSJed Brown 3369566063dSJacob Faibussowitsch PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&A2)); 3379566063dSJacob Faibussowitsch PetscCall(MatDuplicate(S,MAT_COPY_VALUES,&S2)); 3389566063dSJacob Faibussowitsch PetscCall(MatDestroy(&A)); 3399566063dSJacob Faibussowitsch PetscCall(MatDestroy(&S)); 340c4762a1bSJed Brown A = A2; 341c4762a1bSJed Brown S = S2; 342c4762a1bSJed Brown } 343c4762a1bSJed Brown 344c4762a1bSJed Brown if (testsubmat) { 345c4762a1bSJed Brown Mat sA,sS,dA,dS,At,St; 346c4762a1bSJed Brown IS r,c; 347c4762a1bSJed Brown PetscInt rst,ren,cst,cen; 348c4762a1bSJed Brown 3499566063dSJacob Faibussowitsch PetscCall(MatGetOwnershipRange(A,&rst,&ren)); 3509566063dSJacob Faibussowitsch PetscCall(MatGetOwnershipRangeColumn(A,&cst,&cen)); 3519566063dSJacob Faibussowitsch PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A),(ren-rst)/2,rst,1,&r)); 3529566063dSJacob Faibussowitsch PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A),(cen-cst)/2,cst,1,&c)); 3539566063dSJacob Faibussowitsch PetscCall(MatCreateSubMatrix(A,r,c,MAT_INITIAL_MATRIX,&sA)); 3549566063dSJacob Faibussowitsch PetscCall(MatCreateSubMatrix(S,r,c,MAT_INITIAL_MATRIX,&sS)); 3559566063dSJacob Faibussowitsch PetscCall(MatMultAddEqual(sA,sS,10,&flg)); 356c4762a1bSJed Brown if (!flg) { 3579566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error submatrix mult add\n",test)); 358c4762a1bSJed Brown } 3599566063dSJacob Faibussowitsch PetscCall(MatMultTransposeAddEqual(sA,sS,10,&flg)); 360c4762a1bSJed Brown if (!flg) { 3619566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error submatrix mult add (T)\n",test)); 362c4762a1bSJed Brown } 3639566063dSJacob Faibussowitsch PetscCall(MatConvert(sA,MATDENSE,MAT_INITIAL_MATRIX,&dA)); 3649566063dSJacob Faibussowitsch PetscCall(MatConvert(sS,MATDENSE,MAT_INITIAL_MATRIX,&dS)); 3659566063dSJacob Faibussowitsch PetscCall(MatAXPY(dA,-1.0,dS,SAME_NONZERO_PATTERN)); 3669566063dSJacob Faibussowitsch PetscCall(MatNorm(dA,NORM_FROBENIUS,&err)); 367c4762a1bSJed Brown if (err >= tol) { 3689566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error mat submatrix %g\n",test,(double)err)); 369c4762a1bSJed Brown } 3709566063dSJacob Faibussowitsch PetscCall(MatDestroy(&sA)); 3719566063dSJacob Faibussowitsch PetscCall(MatDestroy(&sS)); 3729566063dSJacob Faibussowitsch PetscCall(MatDestroy(&dA)); 3739566063dSJacob Faibussowitsch PetscCall(MatDestroy(&dS)); 3749566063dSJacob Faibussowitsch PetscCall(MatCreateTranspose(A,&At)); 3759566063dSJacob Faibussowitsch PetscCall(MatCreateTranspose(S,&St)); 3769566063dSJacob Faibussowitsch PetscCall(MatCreateSubMatrix(At,c,r,MAT_INITIAL_MATRIX,&sA)); 3779566063dSJacob Faibussowitsch PetscCall(MatCreateSubMatrix(St,c,r,MAT_INITIAL_MATRIX,&sS)); 3789566063dSJacob Faibussowitsch PetscCall(MatMultAddEqual(sA,sS,10,&flg)); 379c4762a1bSJed Brown if (!flg) { 3809566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error submatrix (T) mult add\n",test)); 381c4762a1bSJed Brown } 3829566063dSJacob Faibussowitsch PetscCall(MatMultTransposeAddEqual(sA,sS,10,&flg)); 383c4762a1bSJed Brown if (!flg) { 3849566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error submatrix (T) mult add (T)\n",test)); 385c4762a1bSJed Brown } 3869566063dSJacob Faibussowitsch PetscCall(MatConvert(sA,MATDENSE,MAT_INITIAL_MATRIX,&dA)); 3879566063dSJacob Faibussowitsch PetscCall(MatConvert(sS,MATDENSE,MAT_INITIAL_MATRIX,&dS)); 3889566063dSJacob Faibussowitsch PetscCall(MatAXPY(dA,-1.0,dS,SAME_NONZERO_PATTERN)); 3899566063dSJacob Faibussowitsch PetscCall(MatNorm(dA,NORM_FROBENIUS,&err)); 390c4762a1bSJed Brown if (err >= tol) { 3919566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error mat submatrix (T) %g\n",test,(double)err)); 392c4762a1bSJed Brown } 3939566063dSJacob Faibussowitsch PetscCall(MatDestroy(&sA)); 3949566063dSJacob Faibussowitsch PetscCall(MatDestroy(&sS)); 3959566063dSJacob Faibussowitsch PetscCall(MatDestroy(&dA)); 3969566063dSJacob Faibussowitsch PetscCall(MatDestroy(&dS)); 3979566063dSJacob Faibussowitsch PetscCall(MatDestroy(&At)); 3989566063dSJacob Faibussowitsch PetscCall(MatDestroy(&St)); 3999566063dSJacob Faibussowitsch PetscCall(ISDestroy(&r)); 4009566063dSJacob Faibussowitsch PetscCall(ISDestroy(&c)); 401c4762a1bSJed Brown } 402c4762a1bSJed Brown 403c4762a1bSJed Brown if (testaxpy) { 404c4762a1bSJed Brown Mat tA,tS,dA,dS; 405c4762a1bSJed Brown MatStructure str[3] = { SAME_NONZERO_PATTERN, SUBSET_NONZERO_PATTERN, DIFFERENT_NONZERO_PATTERN }; 406c4762a1bSJed Brown 4079566063dSJacob Faibussowitsch PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&tA)); 408c4762a1bSJed Brown if (testaxpyd && !(test%2)) { 4099566063dSJacob Faibussowitsch PetscCall(PetscObjectReference((PetscObject)tA)); 410c4762a1bSJed Brown tS = tA; 411c4762a1bSJed Brown } else { 4129566063dSJacob Faibussowitsch PetscCall(PetscObjectReference((PetscObject)S)); 413c4762a1bSJed Brown tS = S; 414c4762a1bSJed Brown } 4159566063dSJacob Faibussowitsch PetscCall(MatAXPY(A,0.5,tA,str[test%3])); 4169566063dSJacob Faibussowitsch PetscCall(MatAXPY(S,0.5,tS,str[test%3])); 417c4762a1bSJed Brown /* this will trigger an error the next MatMult or MatMultTranspose call for S */ 4189566063dSJacob Faibussowitsch if (testaxpyerr) PetscCall(MatScale(tA,0)); 4199566063dSJacob Faibussowitsch PetscCall(MatDestroy(&tA)); 4209566063dSJacob Faibussowitsch PetscCall(MatDestroy(&tS)); 4219566063dSJacob Faibussowitsch PetscCall(MatMultAddEqual(A,S,10,&flg)); 422c4762a1bSJed Brown if (!flg) { 4239566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error axpy mult add\n",test)); 424c4762a1bSJed Brown } 4259566063dSJacob Faibussowitsch PetscCall(MatMultTransposeAddEqual(A,S,10,&flg)); 426c4762a1bSJed Brown if (!flg) { 4279566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error axpy mult add (T)\n",test)); 428c4762a1bSJed Brown } 4299566063dSJacob Faibussowitsch PetscCall(MatConvert(A,MATDENSE,MAT_INITIAL_MATRIX,&dA)); 4309566063dSJacob Faibussowitsch PetscCall(MatConvert(S,MATDENSE,MAT_INITIAL_MATRIX,&dS)); 4319566063dSJacob Faibussowitsch PetscCall(MatAXPY(dA,-1.0,dS,SAME_NONZERO_PATTERN)); 4329566063dSJacob Faibussowitsch PetscCall(MatNorm(dA,NORM_FROBENIUS,&err)); 433c4762a1bSJed Brown if (err >= tol) { 4349566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error mat submatrix %g\n",test,(double)err)); 435c4762a1bSJed Brown } 4369566063dSJacob Faibussowitsch PetscCall(MatDestroy(&dA)); 4379566063dSJacob Faibussowitsch PetscCall(MatDestroy(&dS)); 438c4762a1bSJed Brown } 439c4762a1bSJed Brown 440c4762a1bSJed Brown if (testreset && (ntest == 1 || test == ntest-2)) { 441c4762a1bSJed Brown /* reset MATSHELL */ 4429566063dSJacob Faibussowitsch PetscCall(MatAssemblyBegin(S,MAT_FINAL_ASSEMBLY)); 4439566063dSJacob Faibussowitsch PetscCall(MatAssemblyEnd(S,MAT_FINAL_ASSEMBLY)); 444c4762a1bSJed Brown /* reset A */ 4459566063dSJacob Faibussowitsch PetscCall(MatCopy(user->B,A,DIFFERENT_NONZERO_PATTERN)); 446c4762a1bSJed Brown } 447c4762a1bSJed Brown } 448c4762a1bSJed Brown 4499566063dSJacob Faibussowitsch PetscCall(MatDestroy(&A)); 4509566063dSJacob Faibussowitsch PetscCall(MatDestroy(&S)); 4519566063dSJacob Faibussowitsch PetscCall(PetscFree(user)); 4529566063dSJacob Faibussowitsch PetscCall(PetscFinalize()); 453b122ec5aSJacob Faibussowitsch return 0; 454c4762a1bSJed Brown } 455c4762a1bSJed Brown 456c4762a1bSJed Brown /*TEST 457c4762a1bSJed Brown 458c4762a1bSJed Brown testset: 459c4762a1bSJed Brown suffix: rect 460c4762a1bSJed Brown requires: !single 461c4762a1bSJed Brown output_file: output/ex221_1.out 462c4762a1bSJed Brown nsize: {{1 3}} 463c4762a1bSJed Brown args: -loop 3 -keep {{0 1}} -M {{12 19}} -N {{19 12}} -submat {{0 1}} -test_axpy_different {{0 1}} 464c4762a1bSJed Brown 465c4762a1bSJed Brown testset: 466c4762a1bSJed Brown suffix: square 467c4762a1bSJed Brown requires: !single 468c4762a1bSJed Brown output_file: output/ex221_1.out 469c4762a1bSJed Brown nsize: {{1 3}} 470c4762a1bSJed Brown args: -M 21 -N 21 -loop 4 -rows_only {{0 1}} -keep {{0 1}} -submat {{0 1}} -test_axpy_different {{0 1}} 471c4762a1bSJed Brown TEST*/ 472