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 CHKERRQ(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 CHKERRQ(PetscRandomSetInterval(rctx,zero,one)); 24 CHKERRQ(PetscRandomSetFromOptions(rctx)); 25 for (i=0; i<(nr * nc); i++) { 26 CHKERRQ(MatCreateSeqDense(PETSC_COMM_WORLD,arow[i],acol[i],NULL,&subs[i])); 27 } 28 CHKERRQ(MatCreateNest(PETSC_COMM_WORLD,nr,NULL,nc,NULL,subs,&A)); 29 CHKERRQ(MatCreateVecs(A, &x, NULL)); 30 CHKERRQ(MatCreateVecs(A, NULL, &y)); 31 CHKERRQ(VecDuplicate(x, &ATy)); 32 CHKERRQ(VecDuplicate(y, &Ax)); 33 CHKERRQ(MatSetRandom(A,rctx)); 34 CHKERRQ(MatTranspose(A, MAT_INITIAL_MATRIX, &Atranspose)); 35 36 CHKERRQ(MatView(A, PETSC_VIEWER_STDOUT_WORLD)); 37 CHKERRQ(MatNestGetSubMats(A, NULL, NULL, &block)); 38 for (i=0; i<nr; i++) { 39 for (j=0; j<nc; j++) { 40 CHKERRQ(MatView(block[i][j], PETSC_VIEWER_STDOUT_WORLD)); 41 } 42 } 43 44 CHKERRQ(MatView(Atranspose, PETSC_VIEWER_STDOUT_WORLD)); 45 CHKERRQ(MatNestGetSubMats(Atranspose, NULL, NULL, &block)); 46 for (i=0; i<nc; i++) { 47 for (j=0; j<nr; j++) { 48 CHKERRQ(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 CHKERRQ(VecSetRandom(x,rctx)); 55 CHKERRQ(VecSetRandom(y,rctx)); 56 57 CHKERRQ(MatMult(A, x, Ax)); 58 CHKERRQ(VecDot(Ax, y, &dot1)); 59 CHKERRQ(MatMult(Atranspose, y, ATy)); 60 CHKERRQ(VecDot(ATy, x, &dot2)); 61 62 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD, "<Ax, y> = %g\n", (double)PetscRealPart(dot1))); 63 CHKERRQ(PetscPrintf(PETSC_COMM_WORLD, "<x, A^Ty> = %g\n",(double)PetscRealPart(dot2))); 64 } 65 CHKERRQ(VecDestroy(&x)); 66 CHKERRQ(VecDestroy(&y)); 67 CHKERRQ(VecDestroy(&Ax)); 68 69 CHKERRQ(MatCreateSeqDense(PETSC_COMM_WORLD,acol[0]+acol[nr]+acol[2*nr],nk,NULL,&B)); 70 CHKERRQ(MatSetRandom(B,rctx)); 71 CHKERRQ(MatMatMult(A,B,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&C)); 72 CHKERRQ(MatMatMult(A,B,MAT_REUSE_MATRIX,PETSC_DEFAULT,&C)); 73 CHKERRQ(MatMatMultEqual(A,B,C,10,&equal)); 74 PetscCheck(equal,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in C != A*B"); 75 76 CHKERRQ(MatGetSize(A,&M,&N)); 77 CHKERRQ(MatGetLocalSize(A,&m,&n)); 78 for (i=0; i<nk; i++) { 79 CHKERRQ(MatDenseGetColumn(B,i,&valsB)); 80 CHKERRQ(VecCreateMPIWithArray(PETSC_COMM_WORLD,1,n,N,valsB,&x)); 81 CHKERRQ(MatCreateVecs(A,NULL,&Ax)); 82 CHKERRQ(MatMult(A,x,Ax)); 83 CHKERRQ(VecDestroy(&x)); 84 CHKERRQ(MatDenseGetColumn(C,i,&valsC)); 85 CHKERRQ(VecCreateMPIWithArray(PETSC_COMM_WORLD,1,m,M,valsC,&y)); 86 CHKERRQ(VecAXPY(y,-1.0,Ax)); 87 CHKERRQ(VecDestroy(&Ax)); 88 CHKERRQ(VecDestroy(&y)); 89 CHKERRQ(MatDenseRestoreColumn(C,&valsC)); 90 CHKERRQ(MatDenseRestoreColumn(B,&valsB)); 91 } 92 CHKERRQ(MatNorm(C,NORM_INFINITY,&norm)); 93 PetscCheckFalse(norm > PETSC_SMALL,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Error in MatMatMult(): %g",(double)norm); 94 CHKERRQ(MatDestroy(&C)); 95 CHKERRQ(MatDestroy(&B)); 96 97 for (i=0; i<(nr * nc); i++) { 98 CHKERRQ(MatDestroy(&subs[i])); 99 } 100 CHKERRQ(MatDestroy(&A)); 101 CHKERRQ(MatDestroy(&Atranspose)); 102 CHKERRQ(VecDestroy(&ATy)); 103 CHKERRQ(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 CHKERRQ(PetscRandomCreate(PETSC_COMM_WORLD,&rctx)); 118 CHKERRQ(PetscRandomSetInterval(rctx,zero,one)); 119 CHKERRQ(PetscRandomSetFromOptions(rctx)); 120 for (i=0; i<(n * n); i++) { 121 CHKERRQ(MatCreateSeqDense(PETSC_COMM_WORLD,n,n,NULL,&subs[i])); 122 } 123 CHKERRQ(MatCreateNest(PETSC_COMM_WORLD,n,NULL,n,NULL,subs,&A)); 124 CHKERRQ(MatSetRandom(A,rctx)); 125 126 CHKERRQ(MatView(A, PETSC_VIEWER_STDOUT_WORLD)); 127 CHKERRQ(MatNestGetSubMats(A, NULL, NULL, &block)); 128 for (i=0; i<n; i++) { 129 for (j=0; j<n; j++) { 130 CHKERRQ(MatView(block[i][j], PETSC_VIEWER_STDOUT_WORLD)); 131 } 132 } 133 CHKERRQ(MatTranspose(A,MAT_INPLACE_MATRIX,&A)); 134 CHKERRQ(MatView(A, PETSC_VIEWER_STDOUT_WORLD)); 135 CHKERRQ(MatNestGetSubMats(A, NULL, NULL, &block)); 136 for (i=0; i<n; i++) { 137 for (j=0; j<n; j++) { 138 CHKERRQ(MatView(block[i][j], PETSC_VIEWER_STDOUT_WORLD)); 139 } 140 } 141 142 for (i=0; i<(n * n); i++) { 143 CHKERRQ(MatDestroy(&subs[i])); 144 } 145 CHKERRQ(MatDestroy(&A)); 146 CHKERRQ(PetscRandomDestroy(&rctx)); 147 PetscFunctionReturn(0); 148 } 149 150 int main(int argc,char **args) 151 { 152 PetscErrorCode ierr; 153 154 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; 155 CHKERRQ(TestInitialMatrix()); 156 CHKERRQ(TestReuseMatrix()); 157 ierr = PetscFinalize(); 158 return ierr; 159 } 160 161 /*TEST 162 163 test: 164 args: -malloc_dump 165 166 TEST*/ 167