1 static char help[] = "Test MatMatMult(), MatTranspose(), MatTransposeMatMult() for Dense and Elemental matrices.\n\n"; 2 /* 3 Example: 4 mpiexec -n <np> ./ex104 -mat_type elemental 5 */ 6 7 #include <petscmat.h> 8 9 int main(int argc,char **argv) 10 { 11 Mat A,B,C,D; 12 PetscInt i,M=10,N=5,j,nrows,ncols,am,an,rstart,rend; 13 PetscRandom r; 14 PetscBool equal,Aiselemental; 15 PetscReal fill = 1.0; 16 IS isrows,iscols; 17 const PetscInt *rows,*cols; 18 PetscScalar *v,rval; 19 #if defined(PETSC_HAVE_ELEMENTAL) 20 PetscBool Test_MatMatMult=PETSC_TRUE; 21 #else 22 PetscBool Test_MatMatMult=PETSC_FALSE; 23 #endif 24 PetscMPIInt size; 25 26 PetscFunctionBeginUser; 27 PetscCall(PetscInitialize(&argc,&argv,(char*)0,help)); 28 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 29 30 PetscCall(PetscOptionsGetInt(NULL,NULL,"-M",&M,NULL)); 31 PetscCall(PetscOptionsGetInt(NULL,NULL,"-N",&N,NULL)); 32 PetscCall(MatCreate(PETSC_COMM_WORLD,&A)); 33 PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,M,N)); 34 PetscCall(MatSetType(A,MATDENSE)); 35 PetscCall(MatSetFromOptions(A)); 36 PetscCall(MatSetUp(A)); 37 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD,&r)); 38 PetscCall(PetscRandomSetFromOptions(r)); 39 40 /* Set local matrix entries */ 41 PetscCall(MatGetOwnershipIS(A,&isrows,&iscols)); 42 PetscCall(ISGetLocalSize(isrows,&nrows)); 43 PetscCall(ISGetIndices(isrows,&rows)); 44 PetscCall(ISGetLocalSize(iscols,&ncols)); 45 PetscCall(ISGetIndices(iscols,&cols)); 46 PetscCall(PetscMalloc1(nrows*ncols,&v)); 47 for (i=0; i<nrows; i++) { 48 for (j=0; j<ncols; j++) { 49 PetscCall(PetscRandomGetValue(r,&rval)); 50 v[i*ncols+j] = rval; 51 } 52 } 53 PetscCall(MatSetValues(A,nrows,rows,ncols,cols,v,INSERT_VALUES)); 54 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 55 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 56 PetscCall(ISRestoreIndices(isrows,&rows)); 57 PetscCall(ISRestoreIndices(iscols,&cols)); 58 PetscCall(ISDestroy(&isrows)); 59 PetscCall(ISDestroy(&iscols)); 60 PetscCall(PetscRandomDestroy(&r)); 61 62 PetscCall(PetscObjectTypeCompare((PetscObject)A,MATELEMENTAL,&Aiselemental)); 63 64 /* Test MatCreateTranspose() and MatTranspose() */ 65 PetscCall(MatCreateTranspose(A,&C)); 66 PetscCall(MatTranspose(A,MAT_INITIAL_MATRIX,&B)); /* B = A^T */ 67 PetscCall(MatMultEqual(C,B,10,&equal)); 68 PetscCheck(equal,PETSC_COMM_SELF,PETSC_ERR_PLIB,"A^T*x != (x^T*A)^T"); 69 PetscCall(MatDestroy(&B)); 70 71 PetscCall(MatDuplicate(A,MAT_COPY_VALUES,&B)); 72 if (!Aiselemental) { 73 PetscCall(MatTranspose(B,MAT_INPLACE_MATRIX,&B)); 74 PetscCall(MatMultEqual(C,B,10,&equal)); 75 PetscCheck(equal,PETSC_COMM_SELF,PETSC_ERR_PLIB,"C*x != B*x"); 76 } 77 PetscCall(MatDestroy(&B)); 78 79 /* Test B = C*A for matrix type transpose and seqdense */ 80 if (size == 1 && !Aiselemental) { 81 PetscCall(MatMatMult(C,A,MAT_INITIAL_MATRIX,fill,&B)); 82 PetscCall(MatMatMultEqual(C,A,B,10,&equal)); 83 PetscCheck(equal,PETSC_COMM_SELF,PETSC_ERR_PLIB,"B != C*A for matrix type transpose and seqdense"); 84 PetscCall(MatDestroy(&B)); 85 } 86 PetscCall(MatDestroy(&C)); 87 88 /* Test MatMatMult() */ 89 if (Test_MatMatMult) { 90 #if !defined(PETSC_HAVE_ELEMENTAL) 91 PetscCheck(size == 1,PETSC_COMM_SELF,PETSC_ERR_SUP,"This test requires ELEMENTAL"); 92 #endif 93 PetscCall(MatTranspose(A,MAT_INITIAL_MATRIX,&B)); /* B = A^T */ 94 PetscCall(MatMatMult(B,A,MAT_INITIAL_MATRIX,fill,&C)); /* C = B*A = A^T*A */ 95 PetscCall(MatMatMult(B,A,MAT_REUSE_MATRIX,fill,&C)); 96 PetscCall(MatMatMultEqual(B,A,C,10,&equal)); 97 PetscCheck(equal,PETSC_COMM_SELF,PETSC_ERR_PLIB,"B*A*x != C*x"); 98 99 /* Test MatDuplicate for matrix product */ 100 PetscCall(MatDuplicate(C,MAT_COPY_VALUES,&D)); 101 102 PetscCall(MatDestroy(&D)); 103 PetscCall(MatDestroy(&C)); 104 PetscCall(MatDestroy(&B)); 105 } 106 107 /* Test MatTransposeMatMult() */ 108 if (!Aiselemental) { 109 PetscCall(MatTransposeMatMult(A,A,MAT_INITIAL_MATRIX,fill,&D)); /* D = A^T*A */ 110 PetscCall(MatTransposeMatMult(A,A,MAT_REUSE_MATRIX,fill,&D)); 111 PetscCall(MatTransposeMatMultEqual(A,A,D,10,&equal)); 112 PetscCheck(equal,PETSC_COMM_SELF,PETSC_ERR_PLIB,"D*x != A^T*A*x"); 113 114 /* Test MatDuplicate for matrix product */ 115 PetscCall(MatDuplicate(D,MAT_COPY_VALUES,&C)); 116 PetscCall(MatDestroy(&C)); 117 PetscCall(MatDestroy(&D)); 118 119 /* Test D*x = A^T*C*A*x, where C is in AIJ format */ 120 PetscCall(MatGetLocalSize(A,&am,&an)); 121 PetscCall(MatCreate(PETSC_COMM_WORLD,&C)); 122 if (size == 1) { 123 PetscCall(MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,am,am)); 124 } else { 125 PetscCall(MatSetSizes(C,am,am,PETSC_DECIDE,PETSC_DECIDE)); 126 } 127 PetscCall(MatSetFromOptions(C)); 128 PetscCall(MatSetUp(C)); 129 PetscCall(MatGetOwnershipRange(C,&rstart,&rend)); 130 v[0] = 1.0; 131 for (i=rstart; i<rend; i++) { 132 PetscCall(MatSetValues(C,1,&i,1,&i,v,INSERT_VALUES)); 133 } 134 PetscCall(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY)); 135 PetscCall(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY)); 136 137 /* B = C*A, D = A^T*B */ 138 PetscCall(MatMatMult(C,A,MAT_INITIAL_MATRIX,1.0,&B)); 139 PetscCall(MatTransposeMatMult(A,B,MAT_INITIAL_MATRIX,fill,&D)); 140 PetscCall(MatTransposeMatMultEqual(A,B,D,10,&equal)); 141 PetscCheck(equal,PETSC_COMM_SELF,PETSC_ERR_PLIB,"D*x != A^T*B*x"); 142 143 PetscCall(MatDestroy(&D)); 144 PetscCall(MatDestroy(&C)); 145 PetscCall(MatDestroy(&B)); 146 } 147 148 /* Test MatMatTransposeMult() */ 149 if (!Aiselemental) { 150 PetscReal diff, scale; 151 PetscInt am, an, aM, aN; 152 153 PetscCall(MatGetLocalSize(A, &am, &an)); 154 PetscCall(MatGetSize(A, &aM, &aN)); 155 PetscCall(MatCreateDense(PetscObjectComm((PetscObject)A),PETSC_DECIDE, an, aM + 10, aN, NULL, &B)); 156 PetscCall(MatSetRandom(B, NULL)); 157 PetscCall(MatMatTransposeMult(A,B,MAT_INITIAL_MATRIX,fill,&D)); /* D = A*A^T */ 158 159 /* Test MatDuplicate for matrix product */ 160 PetscCall(MatDuplicate(D,MAT_COPY_VALUES,&C)); 161 162 PetscCall(MatMatTransposeMult(A,B,MAT_REUSE_MATRIX,fill,&D)); 163 PetscCall(MatAXPY(C, -1., D, SAME_NONZERO_PATTERN)); 164 165 PetscCall(MatNorm(C, NORM_FROBENIUS, &diff)); 166 PetscCall(MatNorm(D, NORM_FROBENIUS, &scale)); 167 PetscCheck(diff <= PETSC_SMALL * scale,PetscObjectComm((PetscObject)D), PETSC_ERR_PLIB, "MatMatTransposeMult() differs between MAT_INITIAL_MATRIX and MAT_REUSE_MATRIX"); 168 PetscCall(MatDestroy(&C)); 169 170 PetscCall(MatMatTransposeMultEqual(A,B,D,10,&equal)); 171 PetscCheck(equal,PETSC_COMM_SELF,PETSC_ERR_PLIB,"D*x != A^T*A*x"); 172 PetscCall(MatDestroy(&D)); 173 PetscCall(MatDestroy(&B)); 174 175 } 176 177 PetscCall(MatDestroy(&A)); 178 PetscCall(PetscFree(v)); 179 PetscCall(PetscFinalize()); 180 return 0; 181 } 182 183 /*TEST 184 185 test: 186 output_file: output/ex104.out 187 188 test: 189 suffix: 2 190 nsize: 2 191 output_file: output/ex104.out 192 193 test: 194 suffix: 3 195 nsize: 4 196 output_file: output/ex104.out 197 args: -M 23 -N 31 198 199 test: 200 suffix: 4 201 nsize: 4 202 output_file: output/ex104.out 203 args: -M 23 -N 31 -matmattransmult_mpidense_mpidense_via cyclic 204 205 test: 206 suffix: 5 207 nsize: 4 208 output_file: output/ex104.out 209 args: -M 23 -N 31 -matmattransmult_mpidense_mpidense_via allgatherv 210 211 test: 212 suffix: 6 213 args: -mat_type elemental 214 requires: elemental 215 output_file: output/ex104.out 216 217 test: 218 suffix: 7 219 nsize: 2 220 args: -mat_type elemental 221 requires: elemental 222 output_file: output/ex104.out 223 224 TEST*/ 225