1 static char help[] = "Tests the use of MatTranspose_Nest and MatMatMult_Nest_Dense\n"; 2 3 #include <petscmat.h> 4 5 PetscErrorCode TestInitialMatrix(void) 6 { 7 const PetscInt nr = 2,nc = 3,nk = 10; 8 PetscInt n,N,m,M; 9 const PetscInt arow[2*3] = { 2,2,2,3,3,3 }; 10 const PetscInt acol[2*3] = { 3,2,4,3,2,4 }; 11 Mat A,Atranspose,B,C; 12 Mat subs[2*3],**block; 13 Vec x,y,Ax,ATy; 14 PetscInt i,j; 15 PetscScalar dot1,dot2,zero = 0.0,one = 1.0,*valsB,*valsC; 16 PetscReal norm; 17 PetscRandom rctx; 18 PetscBool equal; 19 20 PetscFunctionBegin; 21 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD,&rctx)); 22 /* Force the random numbers to have imaginary part 0 so printed results are the same for --with-scalar-type=real or --with-scalar-type=complex */ 23 PetscCall(PetscRandomSetInterval(rctx,zero,one)); 24 PetscCall(PetscRandomSetFromOptions(rctx)); 25 for (i=0; i<(nr * nc); i++) { 26 PetscCall(MatCreateSeqDense(PETSC_COMM_WORLD,arow[i],acol[i],NULL,&subs[i])); 27 } 28 PetscCall(MatCreateNest(PETSC_COMM_WORLD,nr,NULL,nc,NULL,subs,&A)); 29 PetscCall(MatCreateVecs(A, &x, NULL)); 30 PetscCall(MatCreateVecs(A, NULL, &y)); 31 PetscCall(VecDuplicate(x, &ATy)); 32 PetscCall(VecDuplicate(y, &Ax)); 33 PetscCall(MatSetRandom(A,rctx)); 34 PetscCall(MatTranspose(A, MAT_INITIAL_MATRIX, &Atranspose)); 35 36 PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD)); 37 PetscCall(MatNestGetSubMats(A, NULL, NULL, &block)); 38 for (i=0; i<nr; i++) { 39 for (j=0; j<nc; j++) { 40 PetscCall(MatView(block[i][j], PETSC_VIEWER_STDOUT_WORLD)); 41 } 42 } 43 44 PetscCall(MatView(Atranspose, PETSC_VIEWER_STDOUT_WORLD)); 45 PetscCall(MatNestGetSubMats(Atranspose, NULL, NULL, &block)); 46 for (i=0; i<nc; i++) { 47 for (j=0; j<nr; j++) { 48 PetscCall(MatView(block[i][j], PETSC_VIEWER_STDOUT_WORLD)); 49 } 50 } 51 52 /* Check <Ax, y> = <x, A^Ty> */ 53 for (i=0; i<10; i++) { 54 PetscCall(VecSetRandom(x,rctx)); 55 PetscCall(VecSetRandom(y,rctx)); 56 57 PetscCall(MatMult(A, x, Ax)); 58 PetscCall(VecDot(Ax, y, &dot1)); 59 PetscCall(MatMult(Atranspose, y, ATy)); 60 PetscCall(VecDot(ATy, x, &dot2)); 61 62 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "<Ax, y> = %g\n", (double)PetscRealPart(dot1))); 63 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "<x, A^Ty> = %g\n",(double)PetscRealPart(dot2))); 64 } 65 PetscCall(VecDestroy(&x)); 66 PetscCall(VecDestroy(&y)); 67 PetscCall(VecDestroy(&Ax)); 68 69 PetscCall(MatCreateSeqDense(PETSC_COMM_WORLD,acol[0]+acol[nr]+acol[2*nr],nk,NULL,&B)); 70 PetscCall(MatSetRandom(B,rctx)); 71 PetscCall(MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C)); 72 PetscCall(MatMatMult(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C)); 73 PetscCall(MatMatMultEqual(A,B,C,10,&equal)); 74 PetscCheck(equal,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in C != A*B"); 75 76 PetscCall(MatGetSize(A,&M,&N)); 77 PetscCall(MatGetLocalSize(A,&m,&n)); 78 for (i=0; i<nk; i++) { 79 PetscCall(MatDenseGetColumn(B,i,&valsB)); 80 PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD,1,n,N,valsB,&x)); 81 PetscCall(MatCreateVecs(A,NULL,&Ax)); 82 PetscCall(MatMult(A,x,Ax)); 83 PetscCall(VecDestroy(&x)); 84 PetscCall(MatDenseGetColumn(C,i,&valsC)); 85 PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD,1,m,M,valsC,&y)); 86 PetscCall(VecAXPY(y,-1.0,Ax)); 87 PetscCall(VecDestroy(&Ax)); 88 PetscCall(VecDestroy(&y)); 89 PetscCall(MatDenseRestoreColumn(C,&valsC)); 90 PetscCall(MatDenseRestoreColumn(B,&valsB)); 91 } 92 PetscCall(MatNorm(C,NORM_INFINITY,&norm)); 93 PetscCheckFalse(norm > PETSC_SMALL,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatMatMult(): %g",(double)norm); 94 PetscCall(MatDestroy(&C)); 95 PetscCall(MatDestroy(&B)); 96 97 for (i=0; i<(nr * nc); i++) { 98 PetscCall(MatDestroy(&subs[i])); 99 } 100 PetscCall(MatDestroy(&A)); 101 PetscCall(MatDestroy(&Atranspose)); 102 PetscCall(VecDestroy(&ATy)); 103 PetscCall(PetscRandomDestroy(&rctx)); 104 PetscFunctionReturn(0); 105 } 106 107 PetscErrorCode TestReuseMatrix(void) 108 { 109 const PetscInt n = 2; 110 Mat A; 111 Mat subs[2*2],**block; 112 PetscInt i,j; 113 PetscRandom rctx; 114 PetscScalar zero = 0.0, one = 1.0; 115 116 PetscFunctionBegin; 117 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD,&rctx)); 118 PetscCall(PetscRandomSetInterval(rctx,zero,one)); 119 PetscCall(PetscRandomSetFromOptions(rctx)); 120 for (i=0; i<(n * n); i++) { 121 PetscCall(MatCreateSeqDense(PETSC_COMM_WORLD,n,n,NULL,&subs[i])); 122 } 123 PetscCall(MatCreateNest(PETSC_COMM_WORLD,n,NULL,n,NULL,subs,&A)); 124 PetscCall(MatSetRandom(A,rctx)); 125 126 PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD)); 127 PetscCall(MatNestGetSubMats(A, NULL, NULL, &block)); 128 for (i=0; i<n; i++) { 129 for (j=0; j<n; j++) { 130 PetscCall(MatView(block[i][j], PETSC_VIEWER_STDOUT_WORLD)); 131 } 132 } 133 PetscCall(MatTranspose(A,MAT_INPLACE_MATRIX,&A)); 134 PetscCall(MatView(A, PETSC_VIEWER_STDOUT_WORLD)); 135 PetscCall(MatNestGetSubMats(A, NULL, NULL, &block)); 136 for (i=0; i<n; i++) { 137 for (j=0; j<n; j++) { 138 PetscCall(MatView(block[i][j], PETSC_VIEWER_STDOUT_WORLD)); 139 } 140 } 141 142 for (i=0; i<(n * n); i++) { 143 PetscCall(MatDestroy(&subs[i])); 144 } 145 PetscCall(MatDestroy(&A)); 146 PetscCall(PetscRandomDestroy(&rctx)); 147 PetscFunctionReturn(0); 148 } 149 150 int main(int argc,char **args) 151 { 152 153 PetscCall(PetscInitialize(&argc,&args,(char*)0,help)); 154 PetscCall(TestInitialMatrix()); 155 PetscCall(TestReuseMatrix()); 156 PetscCall(PetscFinalize()); 157 return 0; 158 } 159 160 /*TEST 161 162 test: 163 args: -malloc_dump 164 165 TEST*/ 166