10925cdddSBarry Smith 2e090d566SSatish Balay #include "src/mat/matimpl.h" /*I "petscmat.h" I*/ 30925cdddSBarry Smith 44a2ae208SSatish Balay #undef __FUNCT__ 54a2ae208SSatish Balay #define __FUNCT__ "MatGetColumnVector" 60925cdddSBarry Smith /*@ 782bf6240SBarry Smith MatGetColumnVector - Gets the values from a given column of a matrix. 80925cdddSBarry Smith 972257631SBarry Smith Not Collective 10fee21e36SBarry Smith 1198a79cdbSBarry Smith Input Parameters: 125ed6d96aSBarry Smith + A - the matrix 135ed6d96aSBarry Smith . yy - the vector 145ed6d96aSBarry Smith - c - the column requested (in global numbering) 1598a79cdbSBarry Smith 1615091d37SBarry Smith Level: advanced 1715091d37SBarry Smith 189d006df2SBarry Smith Notes: 199d006df2SBarry Smith Each processor for which this is called gets the values for its rows. 209d006df2SBarry Smith 219d006df2SBarry Smith Since PETSc matrices are usually stored in compressed row format, this routine 229d006df2SBarry Smith will generally be slow. 239d006df2SBarry Smith 243d81755aSBarry Smith The vector must have the same parallel row layout as the matrix. 253d81755aSBarry Smith 2682bf6240SBarry Smith Contributed by: Denis Vanderstraeten 270925cdddSBarry Smith 2882bf6240SBarry Smith .keywords: matrix, column, get 290925cdddSBarry Smith 3015091d37SBarry Smith .seealso: MatGetRow(), MatGetDiagonal() 3115091d37SBarry Smith 320925cdddSBarry Smith @*/ 3382bf6240SBarry Smith int MatGetColumnVector(Mat A,Vec yy,int col) 340925cdddSBarry Smith { 35*b3cc6726SBarry Smith PetscScalar *y,zero = 0.0; 36*b3cc6726SBarry Smith const PetscScalar *v; 37*b3cc6726SBarry Smith int ierr,i,j,nz,N,Rs,Re,rs,re; 38*b3cc6726SBarry Smith const int *idx; 39011b8408SBarry Smith MPI_Comm comm; 400925cdddSBarry Smith 410925cdddSBarry Smith PetscFunctionBegin; 424482741eSBarry Smith PetscValidHeaderSpecific(A,MAT_COOKIE,1); 434482741eSBarry Smith PetscValidHeaderSpecific(yy,VEC_COOKIE,2); 440925cdddSBarry Smith 4529bbc08cSBarry Smith if (col < 0) SETERRQ1(1,"Requested negative column: %d",col); 46011b8408SBarry Smith ierr = MatGetSize(A,PETSC_NULL,&N);CHKERRQ(ierr); 4729bbc08cSBarry Smith if (col >= N) SETERRQ2(1,"Requested column %d larger than number columns in matrix %d",col,N); 480925cdddSBarry Smith 4982bf6240SBarry Smith ierr = MatGetOwnershipRange(A,&Rs,&Re);CHKERRQ(ierr); 50011b8408SBarry Smith 51011b8408SBarry Smith ierr = PetscObjectGetComm((PetscObject)yy,&comm);CHKERRQ(ierr); 5282bf6240SBarry Smith ierr = VecGetOwnershipRange(yy,&rs,&re);CHKERRQ(ierr); 5329bbc08cSBarry Smith 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); 5482bf6240SBarry Smith 5582bf6240SBarry Smith ierr = VecSet(&zero,yy);CHKERRQ(ierr); 5682bf6240SBarry Smith ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 5782bf6240SBarry Smith 5882bf6240SBarry Smith for (i=Rs; i<Re; i++) { 5982bf6240SBarry Smith ierr = MatGetRow(A,i,&nz,&idx,&v);CHKERRQ(ierr); 6082bf6240SBarry Smith if (nz && idx[0] <= col) { 6182bf6240SBarry Smith /* 6282bf6240SBarry Smith Should use faster search here 6382bf6240SBarry Smith */ 6482bf6240SBarry Smith for (j=0; j<nz; j++) { 6582bf6240SBarry Smith if (idx[j] >= col) { 6682bf6240SBarry Smith if (idx[j] == col) y[i-rs] = v[j]; 6782bf6240SBarry Smith break; 680925cdddSBarry Smith } 690925cdddSBarry Smith } 700925cdddSBarry Smith } 7182bf6240SBarry Smith ierr = MatRestoreRow(A,i,&nz,&idx,&v);CHKERRQ(ierr); 720925cdddSBarry Smith } 730925cdddSBarry Smith 7482bf6240SBarry Smith ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 750925cdddSBarry Smith PetscFunctionReturn(0); 760925cdddSBarry Smith } 77