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)); 4708401ef6SPierre 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 79*327415f7SBarry Smith PetscFunctionBeginUser; 809566063dSJacob Faibussowitsch PetscCall(PetscInitialize(&argc,&args,(char*)0,help)); 819566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 829566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 839566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL)); 849566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL)); 859566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetInt(NULL,NULL,"-ml",&m,NULL)); 869566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetInt(NULL,NULL,"-nl",&n,NULL)); 879566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL,NULL,"-square_nc",&nc,NULL)); 889566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL,NULL,"-rows_only",&ronl,NULL)); 899566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL,NULL,"-randomize",&randomize,NULL)); 909566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL,NULL,"-submat",&submat,NULL)); 919566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_zerorows",&testzerorows,NULL)); 929566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_diagscale",&testdiagscale,NULL)); 939566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_getdiag",&testgetdiag,NULL)); 949566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_shift",&testshift,NULL)); 959566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_scale",&testscale,NULL)); 969566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_dup",&testdup,NULL)); 979566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_reset",&testreset,NULL)); 989566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_submat",&testsubmat,NULL)); 999566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_axpy",&testaxpy,NULL)); 1009566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_axpy_different",&testaxpyd,NULL)); 1019566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL,NULL,"-test_axpy_error",&testaxpyerr,NULL)); 1029566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetInt(NULL,NULL,"-loop",&ntest,NULL)); 1039566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetReal(NULL,NULL,"-tol",&tol,NULL)); 1049566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetScalar(NULL,NULL,"-diag",&diag,NULL)); 1059566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetBool(NULL,NULL,"-keep",&keep,NULL)); 106c4762a1bSJed Brown /* This tests square matrices with different row/col layout */ 107c4762a1bSJed Brown if (nc && size > 1) { 108c4762a1bSJed Brown M = PetscMax(PetscMax(N,M),1); 109c4762a1bSJed Brown N = M; 110c4762a1bSJed Brown m = n = 0; 111c4762a1bSJed Brown if (rank == 0) { m = M-1; n = 1; } 112c4762a1bSJed Brown else if (rank == 1) { m = 1; n = N-1; } 113c4762a1bSJed Brown } 1149566063dSJacob Faibussowitsch PetscCall(MatCreateDense(PETSC_COMM_WORLD,m,n,M,N,NULL,&A)); 1159566063dSJacob Faibussowitsch PetscCall(MatGetLocalSize(A,&m,&n)); 1169566063dSJacob Faibussowitsch PetscCall(MatGetSize(A,&M,&N)); 1179566063dSJacob Faibussowitsch PetscCall(MatGetOwnershipRange(A,&s1,NULL)); 118c4762a1bSJed Brown s2 = 1; 119c4762a1bSJed Brown while (s2 < M) s2 *= 10; 1209566063dSJacob Faibussowitsch PetscCall(MatDenseGetArray(A,&data)); 121c4762a1bSJed Brown for (j = 0; j < N; j++) { 122c4762a1bSJed Brown for (i = 0; i < m; i++) { 123c4762a1bSJed Brown data[j*m + i] = s2*j + i + s1 + 1; 124c4762a1bSJed Brown } 125c4762a1bSJed Brown } 1269566063dSJacob Faibussowitsch PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 1279566063dSJacob Faibussowitsch PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 128c4762a1bSJed Brown 129c4762a1bSJed Brown if (submat) { 130c4762a1bSJed Brown Mat A2; 131c4762a1bSJed Brown IS r,c; 132c4762a1bSJed Brown PetscInt rst,ren,cst,cen; 133c4762a1bSJed Brown 1349566063dSJacob Faibussowitsch PetscCall(MatGetOwnershipRange(A,&rst,&ren)); 1359566063dSJacob Faibussowitsch PetscCall(MatGetOwnershipRangeColumn(A,&cst,&cen)); 1369566063dSJacob Faibussowitsch PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A),(ren-rst)/2,rst,1,&r)); 1379566063dSJacob Faibussowitsch PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A),(cen-cst)/2,cst,1,&c)); 1389566063dSJacob Faibussowitsch PetscCall(MatCreateSubMatrix(A,r,c,MAT_INITIAL_MATRIX,&A2)); 1399566063dSJacob Faibussowitsch PetscCall(ISDestroy(&r)); 1409566063dSJacob Faibussowitsch PetscCall(ISDestroy(&c)); 1419566063dSJacob Faibussowitsch PetscCall(MatDestroy(&A)); 142c4762a1bSJed Brown A = A2; 143c4762a1bSJed Brown } 144c4762a1bSJed Brown 1459566063dSJacob Faibussowitsch PetscCall(MatGetSize(A,&M,&N)); 1469566063dSJacob Faibussowitsch PetscCall(MatGetLocalSize(A,&m,&n)); 1479566063dSJacob Faibussowitsch PetscCall(MatHasCongruentLayouts(A,&cong)); 148c4762a1bSJed Brown 1499566063dSJacob Faibussowitsch PetscCall(MatConvert(A,MATAIJ,MAT_INPLACE_MATRIX,&A)); 1509566063dSJacob Faibussowitsch PetscCall(MatSetOption(A,MAT_KEEP_NONZERO_PATTERN,keep)); 1519566063dSJacob Faibussowitsch PetscCall(PetscObjectSetName((PetscObject)A,"initial")); 1529566063dSJacob Faibussowitsch PetscCall(MatViewFromOptions(A,NULL,"-view_mat")); 153c4762a1bSJed Brown 1549566063dSJacob Faibussowitsch PetscCall(PetscNew(&user)); 1559566063dSJacob Faibussowitsch PetscCall(MatCreateShell(PETSC_COMM_WORLD,m,n,M,N,user,&S)); 1569566063dSJacob Faibussowitsch PetscCall(MatShellSetOperation(S,MATOP_MULT,(void (*)(void))MatMult_User)); 1579566063dSJacob Faibussowitsch PetscCall(MatShellSetOperation(S,MATOP_MULT_TRANSPOSE,(void (*)(void))MatMultTranspose_User)); 1581baa6e33SBarry Smith if (cong) PetscCall(MatShellSetOperation(S,MATOP_GET_DIAGONAL,(void (*)(void))MatGetDiagonal_User)); 1599566063dSJacob Faibussowitsch PetscCall(MatShellSetOperation(S,MATOP_COPY,(void (*)(void))MatCopy_User)); 1609566063dSJacob Faibussowitsch PetscCall(MatShellSetOperation(S,MATOP_DESTROY,(void (*)(void))MatDestroy_User)); 1619566063dSJacob Faibussowitsch PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&user->B)); 162c4762a1bSJed Brown 163c4762a1bSJed Brown /* Square and rows only scaling */ 164c4762a1bSJed Brown ronl = cong ? ronl : PETSC_TRUE; 165c4762a1bSJed Brown 166c4762a1bSJed Brown for (test = 0; test < ntest; test++) { 167c4762a1bSJed Brown PetscReal err; 168c4762a1bSJed Brown 1699566063dSJacob Faibussowitsch PetscCall(MatMultAddEqual(A,S,10,&flg)); 170c4762a1bSJed Brown if (!flg) { 1719566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error mult add\n",test)); 172c4762a1bSJed Brown } 1739566063dSJacob Faibussowitsch PetscCall(MatMultTransposeAddEqual(A,S,10,&flg)); 174c4762a1bSJed Brown if (!flg) { 1759566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error mult add (T)\n",test)); 176c4762a1bSJed Brown } 177c4762a1bSJed Brown if (testzerorows) { 178c4762a1bSJed Brown Mat ST,B,C,BT,BTT; 179c4762a1bSJed Brown IS zr; 180c4762a1bSJed Brown Vec x = NULL, b1 = NULL, b2 = NULL; 181c4762a1bSJed Brown PetscInt *idxs = NULL, nr = 0; 182c4762a1bSJed Brown 183c4762a1bSJed Brown if (rank == (test%size)) { 184c4762a1bSJed Brown nr = 1; 1859566063dSJacob Faibussowitsch PetscCall(PetscMalloc1(nr,&idxs)); 186c4762a1bSJed Brown if (test%2) { 187c4762a1bSJed Brown idxs[0] = (2*M - 1 - test/2)%M; 188c4762a1bSJed Brown } else { 189c4762a1bSJed Brown idxs[0] = (test/2)%M; 190c4762a1bSJed Brown } 191c4762a1bSJed Brown idxs[0] = PetscMax(idxs[0],0); 192c4762a1bSJed Brown } 1939566063dSJacob Faibussowitsch PetscCall(ISCreateGeneral(PETSC_COMM_WORLD,nr,idxs,PETSC_OWN_POINTER,&zr)); 1949566063dSJacob Faibussowitsch PetscCall(PetscObjectSetName((PetscObject)zr,"ZR")); 1959566063dSJacob Faibussowitsch PetscCall(ISViewFromOptions(zr,NULL,"-view_is")); 1969566063dSJacob Faibussowitsch PetscCall(MatCreateVecs(A,&x,&b1)); 197c4762a1bSJed Brown if (randomize) { 1989566063dSJacob Faibussowitsch PetscCall(VecSetRandom(x,NULL)); 1999566063dSJacob Faibussowitsch PetscCall(VecSetRandom(b1,NULL)); 200c4762a1bSJed Brown } else { 2019566063dSJacob Faibussowitsch PetscCall(VecSet(x,11.4)); 2029566063dSJacob Faibussowitsch PetscCall(VecSet(b1,-14.2)); 203c4762a1bSJed Brown } 2049566063dSJacob Faibussowitsch PetscCall(VecDuplicate(b1,&b2)); 2059566063dSJacob Faibussowitsch PetscCall(VecCopy(b1,b2)); 2069566063dSJacob Faibussowitsch PetscCall(PetscObjectSetName((PetscObject)b1,"A_B1")); 2079566063dSJacob Faibussowitsch PetscCall(PetscObjectSetName((PetscObject)b2,"A_B2")); 208c4762a1bSJed Brown if (size > 1 && !cong) { /* MATMPIAIJ ZeroRows and ZeroRowsColumns are buggy in this case */ 2099566063dSJacob Faibussowitsch PetscCall(VecDestroy(&b1)); 210c4762a1bSJed Brown } 211c4762a1bSJed Brown if (ronl) { 2129566063dSJacob Faibussowitsch PetscCall(MatZeroRowsIS(A,zr,diag,x,b1)); 2139566063dSJacob Faibussowitsch PetscCall(MatZeroRowsIS(S,zr,diag,x,b2)); 214c4762a1bSJed Brown } else { 2159566063dSJacob Faibussowitsch PetscCall(MatZeroRowsColumnsIS(A,zr,diag,x,b1)); 2169566063dSJacob Faibussowitsch PetscCall(MatZeroRowsColumnsIS(S,zr,diag,x,b2)); 2179566063dSJacob Faibussowitsch PetscCall(ISDestroy(&zr)); 218c4762a1bSJed Brown /* Mix zerorows and zerorowscols */ 219c4762a1bSJed Brown nr = 0; 220c4762a1bSJed Brown idxs = NULL; 221dd400576SPatrick Sanan if (rank == 0) { 222c4762a1bSJed Brown nr = 1; 2239566063dSJacob Faibussowitsch PetscCall(PetscMalloc1(nr,&idxs)); 224c4762a1bSJed Brown if (test%2) { 225c4762a1bSJed Brown idxs[0] = (3*M - 2 - test/2)%M; 226c4762a1bSJed Brown } else { 227c4762a1bSJed Brown idxs[0] = (test/2+1)%M; 228c4762a1bSJed Brown } 229c4762a1bSJed Brown idxs[0] = PetscMax(idxs[0],0); 230c4762a1bSJed Brown } 2319566063dSJacob Faibussowitsch PetscCall(ISCreateGeneral(PETSC_COMM_WORLD,nr,idxs,PETSC_OWN_POINTER,&zr)); 2329566063dSJacob Faibussowitsch PetscCall(PetscObjectSetName((PetscObject)zr,"ZR2")); 2339566063dSJacob Faibussowitsch PetscCall(ISViewFromOptions(zr,NULL,"-view_is")); 2349566063dSJacob Faibussowitsch PetscCall(MatZeroRowsIS(A,zr,diag*2.0+PETSC_SMALL,NULL,NULL)); 2359566063dSJacob Faibussowitsch PetscCall(MatZeroRowsIS(S,zr,diag*2.0+PETSC_SMALL,NULL,NULL)); 236c4762a1bSJed Brown } 2379566063dSJacob Faibussowitsch PetscCall(ISDestroy(&zr)); 238c4762a1bSJed Brown 239c4762a1bSJed Brown if (b1) { 240c4762a1bSJed Brown Vec b; 241c4762a1bSJed Brown 2429566063dSJacob Faibussowitsch PetscCall(VecViewFromOptions(b1,NULL,"-view_b")); 2439566063dSJacob Faibussowitsch PetscCall(VecViewFromOptions(b2,NULL,"-view_b")); 2449566063dSJacob Faibussowitsch PetscCall(VecDuplicate(b1,&b)); 2459566063dSJacob Faibussowitsch PetscCall(VecCopy(b1,b)); 2469566063dSJacob Faibussowitsch PetscCall(VecAXPY(b,-1.0,b2)); 2479566063dSJacob Faibussowitsch PetscCall(VecNorm(b,NORM_INFINITY,&err)); 248c4762a1bSJed Brown if (err >= tol) { 2499566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error b %g\n",test,(double)err)); 250c4762a1bSJed Brown } 2519566063dSJacob Faibussowitsch PetscCall(VecDestroy(&b)); 252c4762a1bSJed Brown } 2539566063dSJacob Faibussowitsch PetscCall(VecDestroy(&b1)); 2549566063dSJacob Faibussowitsch PetscCall(VecDestroy(&b2)); 2559566063dSJacob Faibussowitsch PetscCall(VecDestroy(&x)); 2569566063dSJacob Faibussowitsch PetscCall(MatConvert(S,MATDENSE,MAT_INITIAL_MATRIX,&B)); 257c4762a1bSJed Brown 2589566063dSJacob Faibussowitsch PetscCall(MatCreateTranspose(S,&ST)); 2599566063dSJacob Faibussowitsch PetscCall(MatComputeOperator(ST,MATDENSE,&BT)); 2609566063dSJacob Faibussowitsch PetscCall(MatTranspose(BT,MAT_INITIAL_MATRIX,&BTT)); 2619566063dSJacob Faibussowitsch PetscCall(PetscObjectSetName((PetscObject)B,"S")); 2629566063dSJacob Faibussowitsch PetscCall(PetscObjectSetName((PetscObject)BTT,"STT")); 2639566063dSJacob Faibussowitsch PetscCall(MatConvert(A,MATDENSE,MAT_INITIAL_MATRIX,&C)); 2649566063dSJacob Faibussowitsch PetscCall(PetscObjectSetName((PetscObject)C,"A")); 265c4762a1bSJed Brown 2669566063dSJacob Faibussowitsch PetscCall(MatViewFromOptions(C,NULL,"-view_mat")); 2679566063dSJacob Faibussowitsch PetscCall(MatViewFromOptions(B,NULL,"-view_mat")); 2689566063dSJacob Faibussowitsch PetscCall(MatViewFromOptions(BTT,NULL,"-view_mat")); 269c4762a1bSJed Brown 2709566063dSJacob Faibussowitsch PetscCall(MatAXPY(C,-1.0,B,SAME_NONZERO_PATTERN)); 2719566063dSJacob Faibussowitsch PetscCall(MatNorm(C,NORM_FROBENIUS,&err)); 272c4762a1bSJed Brown if (err >= tol) { 2739566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error mat mult after %s %g\n",test,ronl ? "MatZeroRows" : "MatZeroRowsColumns",(double)err)); 274c4762a1bSJed Brown } 275c4762a1bSJed Brown 2769566063dSJacob Faibussowitsch PetscCall(MatConvert(A,MATDENSE,MAT_REUSE_MATRIX,&C)); 2779566063dSJacob Faibussowitsch PetscCall(MatAXPY(C,-1.0,BTT,SAME_NONZERO_PATTERN)); 2789566063dSJacob Faibussowitsch PetscCall(MatNorm(C,NORM_FROBENIUS,&err)); 279c4762a1bSJed Brown if (err >= tol) { 2809566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error mat mult transpose after %s %g\n",test,ronl ? "MatZeroRows" : "MatZeroRowsColumns",(double)err)); 281c4762a1bSJed Brown } 282c4762a1bSJed Brown 2839566063dSJacob Faibussowitsch PetscCall(MatDestroy(&ST)); 2849566063dSJacob Faibussowitsch PetscCall(MatDestroy(&BTT)); 2859566063dSJacob Faibussowitsch PetscCall(MatDestroy(&BT)); 2869566063dSJacob Faibussowitsch PetscCall(MatDestroy(&B)); 2879566063dSJacob Faibussowitsch PetscCall(MatDestroy(&C)); 288c4762a1bSJed Brown } 289c4762a1bSJed Brown if (testdiagscale) { /* MatDiagonalScale() */ 290c4762a1bSJed Brown Vec vr,vl; 291c4762a1bSJed Brown 2929566063dSJacob Faibussowitsch PetscCall(MatCreateVecs(A,&vr,&vl)); 293c4762a1bSJed Brown if (randomize) { 2949566063dSJacob Faibussowitsch PetscCall(VecSetRandom(vr,NULL)); 2959566063dSJacob Faibussowitsch PetscCall(VecSetRandom(vl,NULL)); 296c4762a1bSJed Brown } else { 2979566063dSJacob Faibussowitsch PetscCall(VecSet(vr,test%2 ? 0.15 : 1.0/0.15)); 2989566063dSJacob Faibussowitsch PetscCall(VecSet(vl,test%2 ? -1.2 : 1.0/-1.2)); 299c4762a1bSJed Brown } 3009566063dSJacob Faibussowitsch PetscCall(MatDiagonalScale(A,vl,vr)); 3019566063dSJacob Faibussowitsch PetscCall(MatDiagonalScale(S,vl,vr)); 3029566063dSJacob Faibussowitsch PetscCall(VecDestroy(&vr)); 3039566063dSJacob Faibussowitsch PetscCall(VecDestroy(&vl)); 304c4762a1bSJed Brown } 305c4762a1bSJed Brown 306c4762a1bSJed Brown if (testscale) { /* MatScale() */ 3079566063dSJacob Faibussowitsch PetscCall(MatScale(A,test%2 ? 1.4 : 1.0/1.4)); 3089566063dSJacob Faibussowitsch PetscCall(MatScale(S,test%2 ? 1.4 : 1.0/1.4)); 309c4762a1bSJed Brown } 310c4762a1bSJed Brown 311c4762a1bSJed Brown if (testshift && cong) { /* MatShift() : MATSHELL shift is broken when row/cols layout are not congruent and left/right scaling have been applied */ 3129566063dSJacob Faibussowitsch PetscCall(MatShift(A,test%2 ? -77.5 : 77.5)); 3139566063dSJacob Faibussowitsch PetscCall(MatShift(S,test%2 ? -77.5 : 77.5)); 314c4762a1bSJed Brown } 315c4762a1bSJed Brown 316c4762a1bSJed Brown if (testgetdiag && cong) { /* MatGetDiagonal() */ 317c4762a1bSJed Brown Vec dA,dS; 318c4762a1bSJed Brown 3199566063dSJacob Faibussowitsch PetscCall(MatCreateVecs(A,&dA,NULL)); 3209566063dSJacob Faibussowitsch PetscCall(MatCreateVecs(S,&dS,NULL)); 3219566063dSJacob Faibussowitsch PetscCall(MatGetDiagonal(A,dA)); 3229566063dSJacob Faibussowitsch PetscCall(MatGetDiagonal(S,dS)); 3239566063dSJacob Faibussowitsch PetscCall(VecAXPY(dA,-1.0,dS)); 3249566063dSJacob Faibussowitsch PetscCall(VecNorm(dA,NORM_INFINITY,&err)); 325c4762a1bSJed Brown if (err >= tol) { 3269566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error diag %g\n",test,(double)err)); 327c4762a1bSJed Brown } 3289566063dSJacob Faibussowitsch PetscCall(VecDestroy(&dA)); 3299566063dSJacob Faibussowitsch PetscCall(VecDestroy(&dS)); 330c4762a1bSJed Brown } 331c4762a1bSJed Brown 332c4762a1bSJed Brown if (testdup && !test) { 333c4762a1bSJed Brown Mat A2, S2; 334c4762a1bSJed Brown 3359566063dSJacob Faibussowitsch PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&A2)); 3369566063dSJacob Faibussowitsch PetscCall(MatDuplicate(S,MAT_COPY_VALUES,&S2)); 3379566063dSJacob Faibussowitsch PetscCall(MatDestroy(&A)); 3389566063dSJacob Faibussowitsch PetscCall(MatDestroy(&S)); 339c4762a1bSJed Brown A = A2; 340c4762a1bSJed Brown S = S2; 341c4762a1bSJed Brown } 342c4762a1bSJed Brown 343c4762a1bSJed Brown if (testsubmat) { 344c4762a1bSJed Brown Mat sA,sS,dA,dS,At,St; 345c4762a1bSJed Brown IS r,c; 346c4762a1bSJed Brown PetscInt rst,ren,cst,cen; 347c4762a1bSJed Brown 3489566063dSJacob Faibussowitsch PetscCall(MatGetOwnershipRange(A,&rst,&ren)); 3499566063dSJacob Faibussowitsch PetscCall(MatGetOwnershipRangeColumn(A,&cst,&cen)); 3509566063dSJacob Faibussowitsch PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A),(ren-rst)/2,rst,1,&r)); 3519566063dSJacob Faibussowitsch PetscCall(ISCreateStride(PetscObjectComm((PetscObject)A),(cen-cst)/2,cst,1,&c)); 3529566063dSJacob Faibussowitsch PetscCall(MatCreateSubMatrix(A,r,c,MAT_INITIAL_MATRIX,&sA)); 3539566063dSJacob Faibussowitsch PetscCall(MatCreateSubMatrix(S,r,c,MAT_INITIAL_MATRIX,&sS)); 3549566063dSJacob Faibussowitsch PetscCall(MatMultAddEqual(sA,sS,10,&flg)); 355c4762a1bSJed Brown if (!flg) { 3569566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error submatrix mult add\n",test)); 357c4762a1bSJed Brown } 3589566063dSJacob Faibussowitsch PetscCall(MatMultTransposeAddEqual(sA,sS,10,&flg)); 359c4762a1bSJed Brown if (!flg) { 3609566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error submatrix mult add (T)\n",test)); 361c4762a1bSJed Brown } 3629566063dSJacob Faibussowitsch PetscCall(MatConvert(sA,MATDENSE,MAT_INITIAL_MATRIX,&dA)); 3639566063dSJacob Faibussowitsch PetscCall(MatConvert(sS,MATDENSE,MAT_INITIAL_MATRIX,&dS)); 3649566063dSJacob Faibussowitsch PetscCall(MatAXPY(dA,-1.0,dS,SAME_NONZERO_PATTERN)); 3659566063dSJacob Faibussowitsch PetscCall(MatNorm(dA,NORM_FROBENIUS,&err)); 366c4762a1bSJed Brown if (err >= tol) { 3679566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error mat submatrix %g\n",test,(double)err)); 368c4762a1bSJed Brown } 3699566063dSJacob Faibussowitsch PetscCall(MatDestroy(&sA)); 3709566063dSJacob Faibussowitsch PetscCall(MatDestroy(&sS)); 3719566063dSJacob Faibussowitsch PetscCall(MatDestroy(&dA)); 3729566063dSJacob Faibussowitsch PetscCall(MatDestroy(&dS)); 3739566063dSJacob Faibussowitsch PetscCall(MatCreateTranspose(A,&At)); 3749566063dSJacob Faibussowitsch PetscCall(MatCreateTranspose(S,&St)); 3759566063dSJacob Faibussowitsch PetscCall(MatCreateSubMatrix(At,c,r,MAT_INITIAL_MATRIX,&sA)); 3769566063dSJacob Faibussowitsch PetscCall(MatCreateSubMatrix(St,c,r,MAT_INITIAL_MATRIX,&sS)); 3779566063dSJacob Faibussowitsch PetscCall(MatMultAddEqual(sA,sS,10,&flg)); 378c4762a1bSJed Brown if (!flg) { 3799566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error submatrix (T) mult add\n",test)); 380c4762a1bSJed Brown } 3819566063dSJacob Faibussowitsch PetscCall(MatMultTransposeAddEqual(sA,sS,10,&flg)); 382c4762a1bSJed Brown if (!flg) { 3839566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error submatrix (T) mult add (T)\n",test)); 384c4762a1bSJed Brown } 3859566063dSJacob Faibussowitsch PetscCall(MatConvert(sA,MATDENSE,MAT_INITIAL_MATRIX,&dA)); 3869566063dSJacob Faibussowitsch PetscCall(MatConvert(sS,MATDENSE,MAT_INITIAL_MATRIX,&dS)); 3879566063dSJacob Faibussowitsch PetscCall(MatAXPY(dA,-1.0,dS,SAME_NONZERO_PATTERN)); 3889566063dSJacob Faibussowitsch PetscCall(MatNorm(dA,NORM_FROBENIUS,&err)); 389c4762a1bSJed Brown if (err >= tol) { 3909566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error mat submatrix (T) %g\n",test,(double)err)); 391c4762a1bSJed Brown } 3929566063dSJacob Faibussowitsch PetscCall(MatDestroy(&sA)); 3939566063dSJacob Faibussowitsch PetscCall(MatDestroy(&sS)); 3949566063dSJacob Faibussowitsch PetscCall(MatDestroy(&dA)); 3959566063dSJacob Faibussowitsch PetscCall(MatDestroy(&dS)); 3969566063dSJacob Faibussowitsch PetscCall(MatDestroy(&At)); 3979566063dSJacob Faibussowitsch PetscCall(MatDestroy(&St)); 3989566063dSJacob Faibussowitsch PetscCall(ISDestroy(&r)); 3999566063dSJacob Faibussowitsch PetscCall(ISDestroy(&c)); 400c4762a1bSJed Brown } 401c4762a1bSJed Brown 402c4762a1bSJed Brown if (testaxpy) { 403c4762a1bSJed Brown Mat tA,tS,dA,dS; 404c4762a1bSJed Brown MatStructure str[3] = { SAME_NONZERO_PATTERN, SUBSET_NONZERO_PATTERN, DIFFERENT_NONZERO_PATTERN }; 405c4762a1bSJed Brown 4069566063dSJacob Faibussowitsch PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&tA)); 407c4762a1bSJed Brown if (testaxpyd && !(test%2)) { 4089566063dSJacob Faibussowitsch PetscCall(PetscObjectReference((PetscObject)tA)); 409c4762a1bSJed Brown tS = tA; 410c4762a1bSJed Brown } else { 4119566063dSJacob Faibussowitsch PetscCall(PetscObjectReference((PetscObject)S)); 412c4762a1bSJed Brown tS = S; 413c4762a1bSJed Brown } 4149566063dSJacob Faibussowitsch PetscCall(MatAXPY(A,0.5,tA,str[test%3])); 4159566063dSJacob Faibussowitsch PetscCall(MatAXPY(S,0.5,tS,str[test%3])); 416c4762a1bSJed Brown /* this will trigger an error the next MatMult or MatMultTranspose call for S */ 4179566063dSJacob Faibussowitsch if (testaxpyerr) PetscCall(MatScale(tA,0)); 4189566063dSJacob Faibussowitsch PetscCall(MatDestroy(&tA)); 4199566063dSJacob Faibussowitsch PetscCall(MatDestroy(&tS)); 4209566063dSJacob Faibussowitsch PetscCall(MatMultAddEqual(A,S,10,&flg)); 421c4762a1bSJed Brown if (!flg) { 4229566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error axpy mult add\n",test)); 423c4762a1bSJed Brown } 4249566063dSJacob Faibussowitsch PetscCall(MatMultTransposeAddEqual(A,S,10,&flg)); 425c4762a1bSJed Brown if (!flg) { 4269566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error axpy mult add (T)\n",test)); 427c4762a1bSJed Brown } 4289566063dSJacob Faibussowitsch PetscCall(MatConvert(A,MATDENSE,MAT_INITIAL_MATRIX,&dA)); 4299566063dSJacob Faibussowitsch PetscCall(MatConvert(S,MATDENSE,MAT_INITIAL_MATRIX,&dS)); 4309566063dSJacob Faibussowitsch PetscCall(MatAXPY(dA,-1.0,dS,SAME_NONZERO_PATTERN)); 4319566063dSJacob Faibussowitsch PetscCall(MatNorm(dA,NORM_FROBENIUS,&err)); 432c4762a1bSJed Brown if (err >= tol) { 4339566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD,"[test %" PetscInt_FMT "] Error mat submatrix %g\n",test,(double)err)); 434c4762a1bSJed Brown } 4359566063dSJacob Faibussowitsch PetscCall(MatDestroy(&dA)); 4369566063dSJacob Faibussowitsch PetscCall(MatDestroy(&dS)); 437c4762a1bSJed Brown } 438c4762a1bSJed Brown 439c4762a1bSJed Brown if (testreset && (ntest == 1 || test == ntest-2)) { 440c4762a1bSJed Brown /* reset MATSHELL */ 4419566063dSJacob Faibussowitsch PetscCall(MatAssemblyBegin(S,MAT_FINAL_ASSEMBLY)); 4429566063dSJacob Faibussowitsch PetscCall(MatAssemblyEnd(S,MAT_FINAL_ASSEMBLY)); 443c4762a1bSJed Brown /* reset A */ 4449566063dSJacob Faibussowitsch PetscCall(MatCopy(user->B,A,DIFFERENT_NONZERO_PATTERN)); 445c4762a1bSJed Brown } 446c4762a1bSJed Brown } 447c4762a1bSJed Brown 4489566063dSJacob Faibussowitsch PetscCall(MatDestroy(&A)); 4499566063dSJacob Faibussowitsch PetscCall(MatDestroy(&S)); 4509566063dSJacob Faibussowitsch PetscCall(PetscFree(user)); 4519566063dSJacob Faibussowitsch PetscCall(PetscFinalize()); 452b122ec5aSJacob Faibussowitsch return 0; 453c4762a1bSJed Brown } 454c4762a1bSJed Brown 455c4762a1bSJed Brown /*TEST 456c4762a1bSJed Brown 457c4762a1bSJed Brown testset: 458c4762a1bSJed Brown suffix: rect 459c4762a1bSJed Brown requires: !single 460c4762a1bSJed Brown output_file: output/ex221_1.out 461c4762a1bSJed Brown nsize: {{1 3}} 462c4762a1bSJed Brown args: -loop 3 -keep {{0 1}} -M {{12 19}} -N {{19 12}} -submat {{0 1}} -test_axpy_different {{0 1}} 463c4762a1bSJed Brown 464c4762a1bSJed Brown testset: 465c4762a1bSJed Brown suffix: square 466c4762a1bSJed Brown requires: !single 467c4762a1bSJed Brown output_file: output/ex221_1.out 468c4762a1bSJed Brown nsize: {{1 3}} 469c4762a1bSJed Brown args: -M 21 -N 21 -loop 4 -rows_only {{0 1}} -keep {{0 1}} -submat {{0 1}} -test_axpy_different {{0 1}} 470c4762a1bSJed Brown TEST*/ 471