Actual source code: ex12.c
petsc-3.12.5 2020-03-29
2: static char help[] = "Reads a PETSc matrix and vector from a file; expands the matrix with the vector\n\n";
4: /*T
5: Concepts: Mat^ordering a matrix - loading a binary matrix and vector;
6: Concepts: Mat^loading a binary matrix and vector;
7: Concepts: Vectors^loading a binary vector;
8: Concepts: PetscLog^preloading executable
9: Processors: 1
10: T*/
14: /*
15: Include "petscmat.h" so that we can use matrices.
16: automatically includes:
17: petscsys.h - base PETSc routines petscvec.h - vectors
18: petscmat.h - matrices
19: petscis.h - index sets petscviewer.h - viewers
20: */
21: #include <petscmat.h>
23: /*
25: Adds a new column and row to the vector (the last) containing the vector
26: */
27: PetscErrorCode PadMatrix(Mat A,Vec v,PetscScalar c,Mat *B)
28: {
29: PetscErrorCode ierr;
30: PetscInt n,i,*cnt,*indices,nc;
31: const PetscInt *aj;
32: const PetscScalar *vv,*aa;
35: MatGetSize(A,&n,NULL);
36: VecGetArrayRead(v,&vv);
37: PetscMalloc1(n,&indices);
38: for (i=0; i<n; i++) indices[i] = i;
40: /* determine number of nonzeros per row in the new matrix */
41: PetscMalloc1(n+1,&cnt);
42: for (i=0; i<n; i++) {
43: MatGetRow(A,i,&nc,NULL,NULL);
44: cnt[i] = nc + (vv[i] != 0.0);
45: MatRestoreRow(A,i,&nc,NULL,NULL);
46: }
47: cnt[n] = 1;
48: for (i=0; i<n; i++) {
49: cnt[n] += (vv[i] != 0.0);
50: }
51: MatCreateSeqAIJ(PETSC_COMM_SELF,n+1,n+1,0,cnt,B);
52: MatSetOption(*B,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);
54: /* copy over the matrix entries from the matrix and then the vector */
55: for (i=0; i<n; i++) {
56: MatGetRow(A,i,&nc,&aj,&aa);
57: MatSetValues(*B,1,&i,nc,aj,aa,INSERT_VALUES);
58: MatRestoreRow(A,i,&nc,&aj,&aa);
59: }
60: MatSetValues(*B,1,&n,n,indices,vv,INSERT_VALUES);
61: MatSetValues(*B,n,indices,1,&n,vv,INSERT_VALUES);
62: MatSetValues(*B,1,&n,1,&n,&c,INSERT_VALUES);
64: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
65: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
66: VecRestoreArrayRead(v,&vv);
67: PetscFree(cnt);
68: PetscFree(indices);
69: return(0);
70: }
73: int main(int argc,char **args)
74: {
75: Mat A,B;
76: PetscViewer fd; /* viewer */
77: char file[PETSC_MAX_PATH_LEN]; /* input file name */
79: PetscBool flg;
80: Vec v;
82: PetscInitialize(&argc,&args,(char*)0,help);if (ierr) return ierr;
83: /*
84: Determine files from which we read the two linear systems
85: (matrix and right-hand-side vector).
86: */
87: PetscOptionsGetString(NULL,NULL,"-f0",file,PETSC_MAX_PATH_LEN,&flg);
88: if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate binary file with the -f0 option");
90: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
92: MatCreate(PETSC_COMM_WORLD,&A);
93: MatSetType(A,MATSEQAIJ);
94: MatLoad(A,fd);
95: VecCreate(PETSC_COMM_WORLD,&v);
96: VecLoad(v,fd);
97: MatView(A,PETSC_VIEWER_STDOUT_SELF);
98: PadMatrix(A,v,3.0,&B);
99: MatView(B,PETSC_VIEWER_STDOUT_SELF);
100: MatDestroy(&B);
101: MatDestroy(&A);
102: VecDestroy(&v);
103: PetscViewerDestroy(&fd);
105: PetscFinalize();
106: return ierr;
107: }
111: /*TEST
113: test:
114: args: -f0 ${wPETSC_DIR}/share/petsc/datafiles/matrices/ns-real-int32-float64
115: requires: double !complex !define(PETSC_USE_64BIT_INDICES)
117: TEST*/