1 2 #include "src/mat/matimpl.h" /*I "petscmat.h" I*/ 3 4 #undef __FUNCT__ 5 #define __FUNCT__ "MatGetColumnVector" 6 /*@ 7 MatGetColumnVector - Gets the values from a given column of a matrix. 8 9 Not Collective 10 11 Input Parameters: 12 + A - the matrix 13 . yy - the vector 14 - c - the column requested (in global numbering) 15 16 Level: advanced 17 18 Notes: 19 Each processor for which this is called gets the values for its rows. 20 21 Since PETSc matrices are usually stored in compressed row format, this routine 22 will generally be slow. 23 24 The vector must have the same parallel row layout as the matrix. 25 26 Contributed by: Denis Vanderstraeten 27 28 .keywords: matrix, column, get 29 30 .seealso: MatGetRow(), MatGetDiagonal() 31 32 @*/ 33 int MatGetColumnVector(Mat A,Vec yy,int col) 34 { 35 PetscScalar *y,zero = 0.0; 36 const PetscScalar *v; 37 int ierr,i,j,nz,N,Rs,Re,rs,re; 38 const int *idx; 39 MPI_Comm comm; 40 41 PetscFunctionBegin; 42 PetscValidHeaderSpecific(A,MAT_COOKIE,1); 43 PetscValidHeaderSpecific(yy,VEC_COOKIE,2); 44 45 if (col < 0) SETERRQ1(1,"Requested negative column: %d",col); 46 ierr = MatGetSize(A,PETSC_NULL,&N);CHKERRQ(ierr); 47 if (col >= N) SETERRQ2(1,"Requested column %d larger than number columns in matrix %d",col,N); 48 49 ierr = MatGetOwnershipRange(A,&Rs,&Re);CHKERRQ(ierr); 50 51 ierr = PetscObjectGetComm((PetscObject)yy,&comm);CHKERRQ(ierr); 52 ierr = VecGetOwnershipRange(yy,&rs,&re);CHKERRQ(ierr); 53 if (Rs != rs || Re != re) SETERRQ4(1,"Matrix %d %d does not have same ownership range (size) as vector %d %d",Rs,Re,rs,re); 54 55 ierr = VecSet(&zero,yy);CHKERRQ(ierr); 56 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 57 58 for (i=Rs; i<Re; i++) { 59 ierr = MatGetRow(A,i,&nz,&idx,&v);CHKERRQ(ierr); 60 if (nz && idx[0] <= col) { 61 /* 62 Should use faster search here 63 */ 64 for (j=0; j<nz; j++) { 65 if (idx[j] >= col) { 66 if (idx[j] == col) y[i-rs] = v[j]; 67 break; 68 } 69 } 70 } 71 ierr = MatRestoreRow(A,i,&nz,&idx,&v);CHKERRQ(ierr); 72 } 73 74 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 75 PetscFunctionReturn(0); 76 } 77