Actual source code: ex12.c
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*/
12: /*
13: Include "petscmat.h" so that we can use matrices.
14: automatically includes:
15: petscsys.h - base PETSc routines petscvec.h - vectors
16: petscmat.h - matrices
17: petscis.h - index sets petscviewer.h - viewers
18: */
19: #include <petscmat.h>
21: /*
23: Adds a new column and row to the vector (the last) containing the vector
24: */
25: PetscErrorCode PadMatrix(Mat A,Vec v,PetscScalar c,Mat *B)
26: {
27: PetscInt n,i,*cnt,*indices,nc;
28: const PetscInt *aj;
29: const PetscScalar *vv,*aa;
31: MatGetSize(A,&n,NULL);
32: VecGetArrayRead(v,&vv);
33: PetscMalloc1(n,&indices);
34: for (i=0; i<n; i++) indices[i] = i;
36: /* determine number of nonzeros per row in the new matrix */
37: PetscMalloc1(n+1,&cnt);
38: for (i=0; i<n; i++) {
39: MatGetRow(A,i,&nc,NULL,NULL);
40: cnt[i] = nc + (vv[i] != 0.0);
41: MatRestoreRow(A,i,&nc,NULL,NULL);
42: }
43: cnt[n] = 1;
44: for (i=0; i<n; i++) {
45: cnt[n] += (vv[i] != 0.0);
46: }
47: MatCreateSeqAIJ(PETSC_COMM_SELF,n+1,n+1,0,cnt,B);
48: MatSetOption(*B,MAT_IGNORE_ZERO_ENTRIES,PETSC_TRUE);
50: /* copy over the matrix entries from the matrix and then the vector */
51: for (i=0; i<n; i++) {
52: MatGetRow(A,i,&nc,&aj,&aa);
53: MatSetValues(*B,1,&i,nc,aj,aa,INSERT_VALUES);
54: MatRestoreRow(A,i,&nc,&aj,&aa);
55: }
56: MatSetValues(*B,1,&n,n,indices,vv,INSERT_VALUES);
57: MatSetValues(*B,n,indices,1,&n,vv,INSERT_VALUES);
58: MatSetValues(*B,1,&n,1,&n,&c,INSERT_VALUES);
60: MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);
61: MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);
62: VecRestoreArrayRead(v,&vv);
63: PetscFree(cnt);
64: PetscFree(indices);
65: return 0;
66: }
68: int main(int argc,char **args)
69: {
70: Mat A,B;
71: PetscViewer fd; /* viewer */
72: char file[PETSC_MAX_PATH_LEN]; /* input file name */
73: PetscBool flg;
74: Vec v;
76: PetscInitialize(&argc,&args,(char*)0,help);
77: /*
78: Determine files from which we read the two linear systems
79: (matrix and right-hand-side vector).
80: */
81: PetscOptionsGetString(NULL,NULL,"-f0",file,sizeof(file),&flg);
84: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,FILE_MODE_READ,&fd);
86: MatCreate(PETSC_COMM_WORLD,&A);
87: MatSetType(A,MATSEQAIJ);
88: MatLoad(A,fd);
89: VecCreate(PETSC_COMM_WORLD,&v);
90: VecLoad(v,fd);
91: MatView(A,PETSC_VIEWER_STDOUT_SELF);
92: PadMatrix(A,v,3.0,&B);
93: MatView(B,PETSC_VIEWER_STDOUT_SELF);
94: MatDestroy(&B);
95: MatDestroy(&A);
96: VecDestroy(&v);
97: PetscViewerDestroy(&fd);
99: PetscFinalize();
100: return 0;
101: }
103: /*TEST
105: test:
106: args: -f0 ${wPETSC_DIR}/share/petsc/datafiles/matrices/ns-real-int32-float64
107: requires: double !complex !defined(PETSC_USE_64BIT_INDICES)
109: TEST*/