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