Actual source code: getcolv.c
petsc-3.11.4 2019-09-28
2: #include <petsc/private/matimpl.h>
4: /*@
5: MatGetColumnVector - Gets the values from a given column of a matrix.
7: Not Collective
9: Input Parameters:
10: + A - the matrix
11: . yy - the vector
12: - col - the column requested (in global numbering)
14: Level: advanced
16: Notes:
17: Each processor for which this is called gets the values for its rows.
19: Since PETSc matrices are usually stored in compressed row format, this routine
20: will generally be slow.
22: The vector must have the same parallel row layout as the matrix.
24: Contributed by: Denis Vanderstraeten
26: .keywords: matrix, column, get
28: .seealso: MatGetRow(), MatGetDiagonal()
30: @*/
31: PetscErrorCode MatGetColumnVector(Mat A,Vec yy,PetscInt col)
32: {
33: PetscScalar *y;
34: const PetscScalar *v;
35: PetscErrorCode ierr;
36: PetscInt i,j,nz,N,Rs,Re,rs,re;
37: const PetscInt *idx;
42: if (col < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Requested negative column: %D",col);
43: MatGetSize(A,NULL,&N);
44: if (col >= N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Requested column %D larger than number columns in matrix %D",col,N);
45: MatGetOwnershipRange(A,&Rs,&Re);
46: VecGetOwnershipRange(yy,&rs,&re);
47: if (Rs != rs || Re != re) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Matrix %D %D does not have same ownership range (size) as vector %D %D",Rs,Re,rs,re);
49: if (A->ops->getcolumnvector) {
50: (*A->ops->getcolumnvector)(A,yy,col);
51: } else {
52: VecSet(yy,0.0);
53: VecGetArray(yy,&y);
55: for (i=Rs; i<Re; i++) {
56: MatGetRow(A,i,&nz,&idx,&v);
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: MatRestoreRow(A,i,&nz,&idx,&v);
69: }
70: VecRestoreArray(yy,&y);
71: }
72: return(0);
73: }
78: /*@
79: MatGetColumnNorms - Gets the norms of each column of a sparse or dense matrix.
81: Input Parameter:
82: + A - the matrix
83: - type - NORM_2, NORM_1 or NORM_INFINITY
85: Output Parameter:
86: . norms - an array as large as the TOTAL number of columns in the matrix
88: Level: intermediate
90: Notes:
91: Each process has ALL the column norms after the call. Because of the way this is computed each process gets all the values,
92: if each process wants only some of the values it should extract the ones it wants from the array.
94: .seealso: NormType, MatNorm()
96: @*/
97: PetscErrorCode MatGetColumnNorms(Mat A,NormType type,PetscReal norms[])
98: {
103: if (A->ops->getcolumnnorms) {
104: (*A->ops->getcolumnnorms)(A,type,norms);
105: } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Not coded for this matrix type");
106: return(0);
107: }