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