1 2 static char help[] = "Tests MatTransposeMatMult() on MatLoad() matrix \n\n"; 3 4 #include <petscmat.h> 5 6 int main(int argc,char **args) 7 { 8 Mat A,C,Bdense,Cdense; 9 PetscErrorCode ierr; 10 PetscViewer fd; /* viewer */ 11 char file[PETSC_MAX_PATH_LEN]; /* input file name */ 12 PetscBool flg,viewmats=PETSC_FALSE; 13 PetscMPIInt rank,size; 14 PetscReal fill=1.0; 15 PetscInt m,n,i,j,BN=10,rstart,rend,*rows,*cols; 16 PetscScalar *Barray,*Carray,rval,*array; 17 Vec x,y; 18 PetscRandom rand; 19 20 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; 21 ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRMPI(ierr); 22 23 /* Determine file from which we read the matrix A */ 24 ierr = PetscOptionsGetString(NULL,NULL,"-f",file,sizeof(file),&flg);CHKERRQ(ierr); 25 if (!flg) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Must indicate binary file with the -f option"); 26 27 /* Load matrix A */ 28 ierr = PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);CHKERRQ(ierr); 29 ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); 30 ierr = MatLoad(A,fd);CHKERRQ(ierr); 31 ierr = PetscViewerDestroy(&fd);CHKERRQ(ierr); 32 33 /* Print (for testing only) */ 34 ierr = PetscOptionsHasName(NULL,NULL, "-view_mats", &viewmats);CHKERRQ(ierr); 35 if (viewmats) { 36 if (!rank) printf("A_aij:\n"); 37 ierr = MatView(A,0);CHKERRQ(ierr); 38 } 39 40 /* Test MatTransposeMatMult_aij_aij() */ 41 ierr = MatTransposeMatMult(A,A,MAT_INITIAL_MATRIX,fill,&C);CHKERRQ(ierr); 42 if (viewmats) { 43 if (!rank) printf("\nC = A_aij^T * A_aij:\n"); 44 ierr = MatView(C,0);CHKERRQ(ierr); 45 } 46 ierr = MatDestroy(&C);CHKERRQ(ierr); 47 ierr = MatGetLocalSize(A,&m,&n);CHKERRQ(ierr); 48 49 /* create a dense matrix Bdense */ 50 ierr = MatCreate(PETSC_COMM_WORLD,&Bdense);CHKERRQ(ierr); 51 ierr = MatSetSizes(Bdense,m,PETSC_DECIDE,PETSC_DECIDE,BN);CHKERRQ(ierr); 52 ierr = MatSetType(Bdense,MATDENSE);CHKERRQ(ierr); 53 ierr = MatSetFromOptions(Bdense);CHKERRQ(ierr); 54 ierr = MatSetUp(Bdense);CHKERRQ(ierr); 55 ierr = MatGetOwnershipRange(Bdense,&rstart,&rend);CHKERRQ(ierr); 56 57 ierr = PetscMalloc3(m,&rows,BN,&cols,m*BN,&array);CHKERRQ(ierr); 58 for (i=0; i<m; i++) rows[i] = rstart + i; 59 ierr = PetscRandomCreate(PETSC_COMM_WORLD,&rand);CHKERRQ(ierr); 60 ierr = PetscRandomSetFromOptions(rand);CHKERRQ(ierr); 61 for (j=0; j<BN; j++) { 62 cols[j] = j; 63 for (i=0; i<m; i++) { 64 ierr = PetscRandomGetValue(rand,&rval);CHKERRQ(ierr); 65 array[m*j+i] = rval; 66 } 67 } 68 ierr = MatSetValues(Bdense,m,rows,BN,cols,array,INSERT_VALUES);CHKERRQ(ierr); 69 ierr = MatAssemblyBegin(Bdense,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 70 ierr = MatAssemblyEnd(Bdense,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 71 ierr = PetscRandomDestroy(&rand);CHKERRQ(ierr); 72 ierr = PetscFree3(rows,cols,array);CHKERRQ(ierr); 73 if (viewmats) { 74 if (!rank) printf("\nBdense:\n"); 75 ierr = MatView(Bdense,0);CHKERRQ(ierr); 76 } 77 78 /* Test MatTransposeMatMult_aij_dense() */ 79 ierr = MatTransposeMatMult(A,Bdense,MAT_INITIAL_MATRIX,fill,&C);CHKERRQ(ierr); 80 ierr = MatTransposeMatMult(A,Bdense,MAT_REUSE_MATRIX,fill,&C);CHKERRQ(ierr); 81 if (viewmats) { 82 if (!rank) printf("\nC=A^T*Bdense:\n"); 83 ierr = MatView(C,0);CHKERRQ(ierr); 84 } 85 86 /* Check accuracy */ 87 ierr = MatCreate(PETSC_COMM_WORLD,&Cdense);CHKERRQ(ierr); 88 ierr = MatSetSizes(Cdense,n,PETSC_DECIDE,PETSC_DECIDE,BN);CHKERRQ(ierr); 89 ierr = MatSetType(Cdense,MATDENSE);CHKERRQ(ierr); 90 ierr = MatSetFromOptions(Cdense);CHKERRQ(ierr); 91 ierr = MatSetUp(Cdense);CHKERRQ(ierr); 92 ierr = MatAssemblyBegin(Cdense,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 93 ierr = MatAssemblyEnd(Cdense,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 94 95 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr); 96 if (size == 1) { 97 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,1,m,NULL,&x);CHKERRQ(ierr); 98 ierr = VecCreateSeqWithArray(PETSC_COMM_SELF,1,n,NULL,&y);CHKERRQ(ierr); 99 } else { 100 ierr = VecCreateMPIWithArray(PETSC_COMM_WORLD,1,m,PETSC_DECIDE,NULL,&x);CHKERRQ(ierr); 101 ierr = VecCreateMPIWithArray(PETSC_COMM_WORLD,1,n,PETSC_DECIDE,NULL,&y);CHKERRQ(ierr); 102 } 103 104 /* Cdense[:,j] = A^T * Bdense[:,j] */ 105 ierr = MatDenseGetArray(Bdense,&Barray);CHKERRQ(ierr); 106 ierr = MatDenseGetArray(Cdense,&Carray);CHKERRQ(ierr); 107 for (j=0; j<BN; j++) { 108 ierr = VecPlaceArray(x,Barray);CHKERRQ(ierr); 109 ierr = VecPlaceArray(y,Carray);CHKERRQ(ierr); 110 111 ierr = MatMultTranspose(A,x,y);CHKERRQ(ierr); 112 113 ierr = VecResetArray(x);CHKERRQ(ierr); 114 ierr = VecResetArray(y);CHKERRQ(ierr); 115 Barray += m; 116 Carray += n; 117 } 118 ierr = MatDenseRestoreArray(Bdense,&Barray);CHKERRQ(ierr); 119 ierr = MatDenseRestoreArray(Cdense,&Carray);CHKERRQ(ierr); 120 if (viewmats) { 121 if (!rank) printf("\nCdense:\n"); 122 ierr = MatView(Cdense,0);CHKERRQ(ierr); 123 } 124 125 ierr = MatEqual(C,Cdense,&flg);CHKERRQ(ierr); 126 if (!flg) { 127 if (!rank) printf(" C != Cdense\n"); 128 } 129 130 /* Free data structures */ 131 ierr = MatDestroy(&A);CHKERRQ(ierr); 132 ierr = MatDestroy(&C);CHKERRQ(ierr); 133 ierr = MatDestroy(&Bdense);CHKERRQ(ierr); 134 ierr = MatDestroy(&Cdense);CHKERRQ(ierr); 135 ierr = VecDestroy(&x);CHKERRQ(ierr); 136 ierr = VecDestroy(&y);CHKERRQ(ierr); 137 ierr = PetscFinalize(); 138 return ierr; 139 } 140 141 /*TEST 142 143 test: 144 requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES) 145 args: -f ${DATAFILESPATH}/matrices/small 146 output_file: output/ex163.out 147 148 test: 149 suffix: 2 150 nsize: 3 151 requires: datafilespath !complex double !define(PETSC_USE_64BIT_INDICES) 152 args: -f ${DATAFILESPATH}/matrices/small 153 output_file: output/ex163.out 154 155 TEST*/ 156