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 followig 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 37 #include <petscmat.h> 38 39 int main(int argc,char **args) 40 { 41 Mat A = NULL; 42 Vec b,u = NULL,u_tmp; 43 char Ain[PETSC_MAX_PATH_LEN],rhs[PETSC_MAX_PATH_LEN],solu[PETSC_MAX_PATH_LEN]; 44 PetscErrorCode ierr; 45 int m,n = 0,nz,dummy; /* these are fscaned so kept as int */ 46 PetscInt i,col,row,shift = 1,sizes[3],nsizes; 47 PetscScalar val; 48 PetscReal res_norm; 49 FILE *Afile,*bfile,*ufile; 50 PetscViewer view; 51 PetscBool flg_A,flg_b,flg_u,flg; 52 PetscMPIInt size; 53 54 ierr = PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr; 55 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr); 56 if (size != 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This is a uniprocessor example only!"); 57 58 /* Read in matrix, rhs and exact solution from ascii files */ 59 ierr = PetscOptionsGetString(NULL,NULL,"-Ain",Ain,PETSC_MAX_PATH_LEN,&flg_A);CHKERRQ(ierr); 60 ierr = PetscOptionsHasName(NULL,NULL,"-noshift",&flg);CHKERRQ(ierr); 61 if (flg) shift = 0; 62 if (flg_A) { 63 ierr = PetscPrintf(PETSC_COMM_SELF,"\n Read matrix in ascii format ...\n");CHKERRQ(ierr); 64 ierr = PetscFOpen(PETSC_COMM_SELF,Ain,"r",&Afile);CHKERRQ(ierr); 65 nsizes = 3; 66 ierr = PetscOptionsGetIntArray(NULL,NULL,"-nosizesinfile",sizes,&nsizes,&flg);CHKERRQ(ierr); 67 if (flg) { 68 if (nsizes != 3) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Must pass in three m,n,nz as arguments for -nosizesinfile"); 69 m = sizes[0]; 70 n = sizes[1]; 71 nz = sizes[2]; 72 } else { 73 if (fscanf(Afile,"%d %d %d\n",&m,&n,&nz) != 3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Badly formatted input file\n"); 74 } 75 ierr = PetscPrintf(PETSC_COMM_SELF,"m: %d, n: %d, nz: %d \n", m,n,nz);CHKERRQ(ierr); 76 if (m != n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ, "Number of rows, cols must be same for this example\n"); 77 ierr = MatCreate(PETSC_COMM_SELF,&A);CHKERRQ(ierr); 78 ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m,n);CHKERRQ(ierr); 79 ierr = MatSetFromOptions(A);CHKERRQ(ierr); 80 ierr = MatSeqAIJSetPreallocation(A,nz/m,NULL);CHKERRQ(ierr); 81 ierr = MatSetOption(A,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);CHKERRQ(ierr); 82 83 for (i=0; i<nz; i++) { 84 if (fscanf(Afile,"%d %d %le\n",&row,&col,(double*)&val) != 3) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Badly formatted input file\n"); 85 row -= shift; col -= shift; /* set index set starts at 0 */ 86 ierr = MatSetValues(A,1,&row,1,&col,&val,INSERT_VALUES);CHKERRQ(ierr); 87 } 88 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 89 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 90 fclose(Afile); 91 } 92 93 ierr = PetscOptionsGetString(NULL,NULL,"-rhs",rhs,PETSC_MAX_PATH_LEN,&flg_b);CHKERRQ(ierr); 94 if (flg_b) { 95 ierr = VecCreate(PETSC_COMM_SELF,&b);CHKERRQ(ierr); 96 ierr = VecSetSizes(b,PETSC_DECIDE,n);CHKERRQ(ierr); 97 ierr = VecSetFromOptions(b);CHKERRQ(ierr); 98 ierr = PetscPrintf(PETSC_COMM_SELF,"\n Read rhs in ascii format ...\n");CHKERRQ(ierr); 99 ierr = PetscFOpen(PETSC_COMM_SELF,rhs,"r",&bfile);CHKERRQ(ierr); 100 for (i=0; i<n; i++) { 101 if (fscanf(bfile,"%d %le\n",&dummy,(double*)&val) != 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Badly formatted input file\n"); 102 ierr = VecSetValues(b,1,&i,&val,INSERT_VALUES);CHKERRQ(ierr); 103 } 104 ierr = VecAssemblyBegin(b);CHKERRQ(ierr); 105 ierr = VecAssemblyEnd(b);CHKERRQ(ierr); 106 fclose(bfile); 107 } 108 109 ierr = PetscOptionsGetString(NULL,NULL,"-solu",solu,PETSC_MAX_PATH_LEN,&flg_u);CHKERRQ(ierr); 110 if (flg_u) { 111 ierr = VecCreate(PETSC_COMM_SELF,&u);CHKERRQ(ierr); 112 ierr = VecSetSizes(u,PETSC_DECIDE,n);CHKERRQ(ierr); 113 ierr = VecSetFromOptions(u);CHKERRQ(ierr); 114 ierr = PetscPrintf(PETSC_COMM_SELF,"\n Read exact solution in ascii format ...\n");CHKERRQ(ierr); 115 ierr = PetscFOpen(PETSC_COMM_SELF,solu,"r",&ufile);CHKERRQ(ierr); 116 for (i=0; i<n; i++) { 117 if (fscanf(ufile,"%d %le\n",&dummy,(double*)&val) != 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Badly formatted input file\n"); 118 ierr = VecSetValues(u,1,&i,&val,INSERT_VALUES);CHKERRQ(ierr); 119 } 120 ierr = VecAssemblyBegin(u);CHKERRQ(ierr); 121 ierr = VecAssemblyEnd(u);CHKERRQ(ierr); 122 fclose(ufile); 123 } 124 125 /* Write matrix, rhs and exact solution in Petsc binary file */ 126 ierr = PetscPrintf(PETSC_COMM_SELF,"\n Write matrix in binary to 'matrix.dat' ...\n");CHKERRQ(ierr); 127 ierr = PetscViewerBinaryOpen(PETSC_COMM_SELF,"matrix.dat",FILE_MODE_WRITE,&view);CHKERRQ(ierr); 128 ierr = MatView(A,view);CHKERRQ(ierr); 129 130 if (flg_b) { /* Write rhs in Petsc binary file */ 131 ierr = PetscPrintf(PETSC_COMM_SELF,"\n Write rhs in binary to 'matrix.dat' ...\n");CHKERRQ(ierr); 132 ierr = VecView(b,view);CHKERRQ(ierr); 133 } 134 if (flg_u) { 135 ierr = PetscPrintf(PETSC_COMM_SELF,"\n Write exact solution in binary to 'matrix.dat' ...\n");CHKERRQ(ierr); 136 ierr = VecView(u,view);CHKERRQ(ierr); 137 } 138 139 /* Check accuracy of the data */ 140 if (flg_A & flg_b & flg_u) { 141 ierr = VecDuplicate(u,&u_tmp);CHKERRQ(ierr); 142 ierr = MatMult(A,u,u_tmp);CHKERRQ(ierr); 143 ierr = VecAXPY(u_tmp,-1.0,b);CHKERRQ(ierr); 144 ierr = VecNorm(u_tmp,NORM_2,&res_norm);CHKERRQ(ierr); 145 ierr = PetscPrintf(PETSC_COMM_SELF,"\n Accuracy of the reading data: | b - A*u |_2 : %g \n",res_norm);CHKERRQ(ierr); 146 ierr = VecDestroy(&u_tmp);CHKERRQ(ierr); 147 } 148 149 ierr = MatDestroy(&A);CHKERRQ(ierr); 150 if (flg_b) {ierr = VecDestroy(&b);CHKERRQ(ierr);} 151 if (flg_u) {ierr = VecDestroy(&u);CHKERRQ(ierr);} 152 ierr = PetscViewerDestroy(&view);CHKERRQ(ierr); 153 ierr = PetscFinalize(); 154 return ierr; 155 } 156 157 /*TEST 158 159 build: 160 requires: !define(PETSC_USE_64BIT_INDICES) double !complex datafilespath 161 162 test: 163 args: -Ain ${DATAFILESPATH}/matrices/indefinite/afiro_A.dat 164 165 TEST*/ 166