Actual source code: ex4.c
1: /*$Id: ex1.c,v 1.28 2001/04/10 19:35:58 bsmith Exp $*/
3: static char help[] = "Reads U and V matrices from a file and performs y = V*U'*x.\n\
4: -f <input_file> : file to load \n\n";
6: /*
7: Include "petscmat.h" so that we can use matrices.
8: automatically includes:
9: petsc.h - base PETSc routines petscvec.h - vectors
10: petscsys.h - system routines petscmat.h - matrices
11: petscis.h - index sets petscviewer.h - viewers
12: */
13: #include petscmat.h
14: extern int LowRankUpdate(Mat,Mat,Vec,Vec,Vec,Vec,int);
19: int main(int argc,char **args)
20: {
21: Mat U,V; /* matrix */
22: PetscViewer fd; /* viewer */
23: char file[128]; /* input file name */
24: int ierr;
25: PetscTruth flg;
26: Vec x,y,work1,work2;
27: int i,N,n,M,m;
28: PetscScalar *xx;
30: PetscInitialize(&argc,&args,(char *)0,help);
32: /*
33: Determine file from which we read the matrix
35: */
36: PetscOptionsGetString(PETSC_NULL,"-f",file,127,&flg);
37: if (!flg) SETERRQ(1,"Must indicate binary file with the -f option");
40: /*
41: Open binary file. Note that we use PETSC_FILE_RDONLY to indicate
42: reading from this file.
43: */
44: PetscViewerBinaryOpen(PETSC_COMM_WORLD,file,PETSC_FILE_RDONLY,&fd);
46: /*
47: Load the matrix; then destroy the viewer.
48: Note both U and V are stored as tall skinny matrices
49: */
50: MatLoad(fd,MATMPIDENSE,&U);
51: MatLoad(fd,MATMPIDENSE,&V);
52: PetscViewerDestroy(fd);
54: MatGetLocalSize(U,&N,&n);
55: MatGetLocalSize(V,&M,&m);
56: if (N != M) SETERRQ2(1,"U and V matrices must have same number of local rows %d %d",N,M);
57: if (n != m) SETERRQ2(1,"U and V matrices must have same number of local columns %d %d",n,m);
59: VecCreateMPI(PETSC_COMM_WORLD,N,PETSC_DETERMINE,&x);
60: VecDuplicate(x,&y);
62: MatGetSize(U,0,&n);
63: VecCreateSeq(PETSC_COMM_SELF,n,&work1);
64: VecDuplicate(work1,&work2);
66: /* put some initial values into x for testing */
67: VecGetArray(x,&xx);
68: for (i=0; i<N; i++) xx[i] = i;
69: VecRestoreArray(x,&xx);
70: LowRankUpdate(U,V,x,y,work1,work2,n);
71: VecView(x,PETSC_VIEWER_STDOUT_WORLD);
72: VecView(y,PETSC_VIEWER_STDOUT_WORLD);
74: /*
75: Free work space. All PETSc objects should be destroyed when they
76: are no longer needed.
77: */
78: MatDestroy(U);
79: MatDestroy(V);
80: VecDestroy(x);
81: VecDestroy(y);
82: VecDestroy(work1);
83: VecDestroy(work2);
85: PetscFinalize();
86: return 0;
87: }
89: #include src/mat/impls/dense/mpi/mpidense.h
93: /*
94: Computes y = V*U'*x where U and V are N by n (N >> n).
96: U and V are stored as PETSc MPIDENSE (parallel) dense matrices with their rows partitioned the
97: same way as x and y are partitioned
98: */
99: int LowRankUpdate(Mat U,Mat V,Vec x,Vec y,Vec work1,Vec work2,int nwork)
100: {
101: Mat Ulocal = ((Mat_MPIDense*)U->data)->A;
102: Mat Vlocal = ((Mat_MPIDense*)V->data)->A;
103: int ierr,Nsave = x->N;
104: PetscScalar *w1,*w2;
108: /* First multiply the local part of U with the local part of x */
109: x->N = x->n; /* this tricks the silly error checking in MatMultTranspose();*/
110: MatMultTranspose(Ulocal,x,work1);/* note in this call x is treated as a sequential vector */
111: x->N = Nsave;
113: /* Form the sum of all the local multiplies : this is work2 = U'*x = sum_{all processors} work1 */
114: VecGetArray(work1,&w1);
115: VecGetArray(work2,&w2);
116: MPI_Allreduce(w1,w2,nwork,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD);
117: VecRestoreArray(work1,&w1);
118: VecRestoreArray(work2,&w2);
120: /* multiply y = V*work2 */
121: y->N = y->n; /* this tricks the silly error checking in MatMult() */
122: MatMult(Vlocal,work2,y);/* note in this call y is treated as a sequential vector */
123: y->N = Nsave;
125: return(0);
126: }