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