1 2 static char help[] = "Reads in a matrix in ASCII MATLAB format (I,J,A), read in vectors rhs and exact_solu in ASCII format.\n\ 3 Writes them using the PETSc sparse format.\n\ 4 Note: I and J start at 1, not 0, use -noshift if indices in file start with zero!\n\ 5 Input parameters are:\n\ 6 -Ain <filename> : input matrix in ascii format\n\ 7 -rhs <filename> : input rhs in ascii format\n\ 8 -solu <filename> : input true solution in ascii format\n\\n"; 9 10 /* 11 Example: ./ex78 -Ain Ain -rhs rhs -solu solu -noshift -mat_view 12 with the datafiles in the following format: 13 Ain (I and J start at 0): 14 ------------------------ 15 3 3 6 16 0 0 1.0 17 0 1 2.0 18 1 0 3.0 19 1 1 4.0 20 1 2 5.0 21 2 2 6.0 22 23 rhs 24 --- 25 0 3.0 26 1 12.0 27 2 6.0 28 29 solu 30 ---- 31 0 1.0 32 0 1.0 33 0 1.0 34 */ 35 36 #include <petscmat.h> 37 38 int main(int argc,char **args) 39 { 40 Mat A = NULL; 41 Vec b,u = NULL,u_tmp; 42 char Ain[PETSC_MAX_PATH_LEN],rhs[PETSC_MAX_PATH_LEN],solu[PETSC_MAX_PATH_LEN]; 43 int m,n = 0,nz,dummy; /* these are fscaned so kept as int */ 44 PetscInt i,col,row,shift = 1,sizes[3],nsizes; 45 PetscScalar val; 46 PetscReal res_norm; 47 FILE *Afile,*bfile,*ufile; 48 PetscViewer view; 49 PetscBool flg_A,flg_b,flg_u,flg; 50 PetscMPIInt size; 51 52 PetscFunctionBeginUser; 53 PetscCall(PetscInitialize(&argc,&args,(char*)0,help)); 54 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 55 PetscCheck(size == 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This is a uniprocessor example only!"); 56 57 /* Read in matrix, rhs and exact solution from ascii files */ 58 PetscCall(PetscOptionsGetString(NULL,NULL,"-Ain",Ain,sizeof(Ain),&flg_A)); 59 PetscCall(PetscOptionsHasName(NULL,NULL,"-noshift",&flg)); 60 if (flg) shift = 0; 61 if (flg_A) { 62 PetscCall(PetscPrintf(PETSC_COMM_SELF,"\n Read matrix in ascii format ...\n")); 63 PetscCall(PetscFOpen(PETSC_COMM_SELF,Ain,"r",&Afile)); 64 nsizes = 3; 65 PetscCall(PetscOptionsGetIntArray(NULL,NULL,"-nosizesinfile",sizes,&nsizes,&flg)); 66 if (flg) { 67 PetscCheck(nsizes == 3,PETSC_COMM_WORLD,PETSC_ERR_USER,"Must pass in three m,n,nz as arguments for -nosizesinfile"); 68 m = sizes[0]; 69 n = sizes[1]; 70 nz = sizes[2]; 71 } else { 72 PetscCheck(fscanf(Afile,"%d %d %d\n",&m,&n,&nz) == 3,PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Badly formatted input file"); 73 } 74 PetscCall(PetscPrintf(PETSC_COMM_SELF,"m: %d, n: %d, nz: %d \n", m,n,nz)); 75 PetscCheck(m == n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "Number of rows, cols must be same for this example"); 76 PetscCall(MatCreate(PETSC_COMM_SELF,&A)); 77 PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m,n)); 78 PetscCall(MatSetFromOptions(A)); 79 PetscCall(MatSeqAIJSetPreallocation(A,nz/m,NULL)); 80 PetscCall(MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE)); 81 82 for (i=0; i<nz; i++) { 83 PetscCheck(fscanf(Afile,"%d %d %le\n",&row,&col,(double*)&val) == 3,PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Badly formatted input file"); 84 row -= shift; col -= shift; /* set index set starts at 0 */ 85 PetscCall(MatSetValues(A,1,&row,1,&col,&val,INSERT_VALUES)); 86 } 87 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 88 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 89 fclose(Afile); 90 } 91 92 PetscCall(PetscOptionsGetString(NULL,NULL,"-rhs",rhs,sizeof(rhs),&flg_b)); 93 if (flg_b) { 94 PetscCall(VecCreate(PETSC_COMM_SELF,&b)); 95 PetscCall(VecSetSizes(b,PETSC_DECIDE,n)); 96 PetscCall(VecSetFromOptions(b)); 97 PetscCall(PetscPrintf(PETSC_COMM_SELF,"\n Read rhs in ascii format ...\n")); 98 PetscCall(PetscFOpen(PETSC_COMM_SELF,rhs,"r",&bfile)); 99 for (i=0; i<n; i++) { 100 PetscCheck(fscanf(bfile,"%d %le\n",&dummy,(double*)&val) == 2,PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Badly formatted input file"); 101 PetscCall(VecSetValues(b,1,&i,&val,INSERT_VALUES)); 102 } 103 PetscCall(VecAssemblyBegin(b)); 104 PetscCall(VecAssemblyEnd(b)); 105 fclose(bfile); 106 } 107 108 PetscCall(PetscOptionsGetString(NULL,NULL,"-solu",solu,sizeof(solu),&flg_u)); 109 if (flg_u) { 110 PetscCall(VecCreate(PETSC_COMM_SELF,&u)); 111 PetscCall(VecSetSizes(u,PETSC_DECIDE,n)); 112 PetscCall(VecSetFromOptions(u)); 113 PetscCall(PetscPrintf(PETSC_COMM_SELF,"\n Read exact solution in ascii format ...\n")); 114 PetscCall(PetscFOpen(PETSC_COMM_SELF,solu,"r",&ufile)); 115 for (i=0; i<n; i++) { 116 PetscCheck(fscanf(ufile,"%d %le\n",&dummy,(double*)&val) == 2,PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Badly formatted input file"); 117 PetscCall(VecSetValues(u,1,&i,&val,INSERT_VALUES)); 118 } 119 PetscCall(VecAssemblyBegin(u)); 120 PetscCall(VecAssemblyEnd(u)); 121 fclose(ufile); 122 } 123 124 /* Write matrix, rhs and exact solution in Petsc binary file */ 125 PetscCall(PetscPrintf(PETSC_COMM_SELF,"\n Write matrix in binary to 'matrix.dat' ...\n")); 126 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_SELF,"matrix.dat",FILE_MODE_WRITE,&view)); 127 PetscCall(MatView(A,view)); 128 129 if (flg_b) { /* Write rhs in Petsc binary file */ 130 PetscCall(PetscPrintf(PETSC_COMM_SELF,"\n Write rhs in binary to 'matrix.dat' ...\n")); 131 PetscCall(VecView(b,view)); 132 } 133 if (flg_u) { 134 PetscCall(PetscPrintf(PETSC_COMM_SELF,"\n Write exact solution in binary to 'matrix.dat' ...\n")); 135 PetscCall(VecView(u,view)); 136 } 137 138 /* Check accuracy of the data */ 139 if (flg_A & flg_b & flg_u) { 140 PetscCall(VecDuplicate(u,&u_tmp)); 141 PetscCall(MatMult(A,u,u_tmp)); 142 PetscCall(VecAXPY(u_tmp,-1.0,b)); 143 PetscCall(VecNorm(u_tmp,NORM_2,&res_norm)); 144 PetscCall(PetscPrintf(PETSC_COMM_SELF,"\n Accuracy of the reading data: | b - A*u |_2 : %g \n",res_norm)); 145 PetscCall(VecDestroy(&u_tmp)); 146 } 147 148 PetscCall(MatDestroy(&A)); 149 if (flg_b) PetscCall(VecDestroy(&b)); 150 if (flg_u) PetscCall(VecDestroy(&u)); 151 PetscCall(PetscViewerDestroy(&view)); 152 PetscCall(PetscFinalize()); 153 return 0; 154 } 155 156 /*TEST 157 158 build: 159 requires: !defined(PETSC_USE_64BIT_INDICES) double !complex 160 161 test: 162 requires: datafilespath 163 args: -Ain ${DATAFILESPATH}/matrices/indefinite/afiro_A.dat 164 165 TEST*/ 166