1 /*$Id: getcolv.c,v 1.14 2000/05/04 03:12:33 bsmith Exp balay $*/ 2 3 #include "src/mat/matimpl.h" /*I "petscmat.h" I*/ 4 5 #undef __FUNC__ 6 #define __FUNC__ /*<a name=""></a>*/"MatGetColumnVector" 7 /*@ 8 MatGetColumnVector - Gets the values from a given column of a matrix. 9 10 Not Collective 11 12 Input Parameters: 13 + X - the matrix 14 . v - the vector 15 - c - the column requested 16 17 Level: advanced 18 19 Contributed by: Denis Vanderstraeten 20 21 .keywords: matrix, column, get 22 23 .seealso: MatGetRow(), MatGetDiagonal() 24 25 @*/ 26 int MatGetColumnVector(Mat A,Vec yy,int col) 27 { 28 Scalar *y,*v,zero = 0.0; 29 int ierr,i,j,nz,*idx,N,Rs,Re,rs,re,size,mlocal; 30 MPI_Comm comm; 31 32 PetscFunctionBegin; 33 PetscValidHeaderSpecific(A,MAT_COOKIE); 34 PetscValidHeaderSpecific(yy,VEC_COOKIE); 35 36 if (col < 0) SETERRQ1(1,1,"Requested negative column: %d",col); 37 ierr = MatGetSize(A,PETSC_NULL,&N);CHKERRQ(ierr); 38 if (col >= N) SETERRQ2(1,1,"Requested column %d larger than number columns in matrix %d",col,N); 39 40 ierr = MatGetOwnershipRange(A,&Rs,&Re);CHKERRQ(ierr); 41 42 ierr = PetscObjectGetComm((PetscObject)yy,&comm);CHKERRQ(ierr); 43 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 44 if (size > 1) { 45 ierr = VecGetOwnershipRange(yy,&rs,&re);CHKERRQ(ierr); 46 if (Rs != rs || Re != re) SETERRQ4(1,1,"Matrix %d %d does not have same ownership range as parallel vector %d %d",Rs,Re,rs,re); 47 } else { 48 ierr = VecGetSize(yy,&mlocal);CHKERRQ(ierr); 49 if (mlocal != Re - Rs) SETERRQ2(1,1,"Matrix %d does not have same ownership size as vector %d",Re-Rs,mlocal); 50 } 51 52 ierr = VecSet(&zero,yy);CHKERRQ(ierr); 53 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 54 55 for (i=Rs; i<Re; i++) { 56 ierr = MatGetRow(A,i,&nz,&idx,&v);CHKERRQ(ierr); 57 if (nz && idx[0] <= col) { 58 /* 59 Should use faster search here 60 */ 61 for (j=0; j<nz; j++) { 62 if (idx[j] >= col) { 63 if (idx[j] == col) y[i-rs] = v[j]; 64 break; 65 } 66 } 67 } 68 ierr = MatRestoreRow(A,i,&nz,&idx,&v);CHKERRQ(ierr); 69 } 70 71 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 72 PetscFunctionReturn(0); 73 } 74