1 2 #include "private/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 PetscErrorCode MatGetColumnVector(Mat A,Vec yy,PetscInt col) 34 { 35 PetscScalar *y; 36 const PetscScalar *v; 37 PetscErrorCode ierr; 38 PetscInt i,j,nz,N,Rs,Re,rs,re; 39 const PetscInt *idx; 40 41 PetscFunctionBegin; 42 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 43 PetscValidHeaderSpecific(yy,VEC_CLASSID,2); 44 if (col < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Requested negative column: %D",col); 45 ierr = MatGetSize(A,PETSC_NULL,&N);CHKERRQ(ierr); 46 if (col >= N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Requested column %D larger than number columns in matrix %D",col,N); 47 ierr = MatGetOwnershipRange(A,&Rs,&Re);CHKERRQ(ierr); 48 ierr = VecGetOwnershipRange(yy,&rs,&re);CHKERRQ(ierr); 49 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); 50 51 if (A->ops->getcolumnvector) { 52 ierr = (*A->ops->getcolumnvector)(A,yy,col);CHKERRQ(ierr); 53 } else { 54 ierr = VecSet(yy,0.0);CHKERRQ(ierr); 55 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 56 57 for (i=Rs; i<Re; i++) { 58 ierr = MatGetRow(A,i,&nz,&idx,&v);CHKERRQ(ierr); 59 if (nz && idx[0] <= col) { 60 /* 61 Should use faster search here 62 */ 63 for (j=0; j<nz; j++) { 64 if (idx[j] >= col) { 65 if (idx[j] == col) y[i-rs] = v[j]; 66 break; 67 } 68 } 69 } 70 ierr = MatRestoreRow(A,i,&nz,&idx,&v);CHKERRQ(ierr); 71 } 72 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 73 } 74 PetscFunctionReturn(0); 75 } 76 77 #include "../src/mat/impls/aij/seq/aij.h" 78 79 #undef __FUNCT__ 80 #define __FUNCT__ "MatGetColumnNorms_SeqAIJ" 81 PetscErrorCode MatGetColumnNorms_SeqAIJ(Mat A,NormType type,PetscReal *norms) 82 { 83 PetscErrorCode ierr; 84 PetscInt i,m,n; 85 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 86 87 PetscFunctionBegin; 88 ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr); 89 ierr = PetscMemzero(norms,n*sizeof(PetscReal));CHKERRQ(ierr); 90 if (type == NORM_2) { 91 for (i=0; i<aij->i[m]; i++) { 92 norms[aij->j[i]] += PetscAbsScalar(aij->a[i]*aij->a[i]); 93 } 94 } else if (type == NORM_1) { 95 for (i=0; i<aij->i[m]; i++) { 96 norms[aij->j[i]] += PetscAbsScalar(aij->a[i]); 97 } 98 } else if (type == NORM_INFINITY) { 99 for (i=0; i<aij->i[m]; i++) { 100 norms[aij->j[i]] = PetscMax(PetscAbsScalar(aij->a[i]),norms[aij->j[i]]); 101 } 102 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unknown NormType"); 103 104 if (type == NORM_2) { 105 for (i=0; i<n; i++) norms[i] = sqrt(norms[i]); 106 } 107 PetscFunctionReturn(0); 108 } 109 110 #include "../src/mat/impls/aij/mpi/mpiaij.h" 111 112 #undef __FUNCT__ 113 #define __FUNCT__ "MatGetColumnNorms_MPIAIJ" 114 PetscErrorCode MatGetColumnNorms_MPIAIJ(Mat A,NormType type,PetscReal *norms) 115 { 116 PetscErrorCode ierr; 117 Mat_MPIAIJ *aij = (Mat_MPIAIJ*)A->data; 118 PetscInt i,n,*garray = aij->garray; 119 Mat_SeqAIJ *a_aij = (Mat_SeqAIJ*) aij->A->data; 120 Mat_SeqAIJ *b_aij = (Mat_SeqAIJ*) aij->B->data; 121 PetscReal *work; 122 123 PetscFunctionBegin; 124 ierr = MatGetSize(A,PETSC_NULL,&n);CHKERRQ(ierr); 125 ierr = PetscMalloc(n*sizeof(PetscReal),&work);CHKERRQ(ierr); 126 ierr = PetscMemzero(work,n*sizeof(PetscReal));CHKERRQ(ierr); 127 if (type == NORM_2) { 128 for (i=0; i<a_aij->i[aij->A->rmap->n]; i++) { 129 work[A->cmap->rstart + a_aij->j[i]] += PetscAbsScalar(a_aij->a[i]*a_aij->a[i]); 130 } 131 for (i=0; i<b_aij->i[aij->B->rmap->n]; i++) { 132 work[garray[b_aij->j[i]]] += PetscAbsScalar(b_aij->a[i]*b_aij->a[i]); 133 } 134 } else if (type == NORM_1) { 135 for (i=0; i<a_aij->i[aij->A->rmap->n]; i++) { 136 work[A->cmap->rstart + a_aij->j[i]] += PetscAbsScalar(a_aij->a[i]); 137 } 138 for (i=0; i<b_aij->i[aij->B->rmap->n]; i++) { 139 work[garray[b_aij->j[i]]] += PetscAbsScalar(b_aij->a[i]); 140 } 141 } else if (type == NORM_INFINITY) { 142 for (i=0; i<a_aij->i[aij->A->rmap->n]; i++) { 143 work[A->cmap->rstart + a_aij->j[i]] = PetscMax(PetscAbsScalar(a_aij->a[i]), work[A->cmap->rstart + a_aij->j[i]]); 144 } 145 for (i=0; i<b_aij->i[aij->B->rmap->n]; i++) { 146 work[garray[b_aij->j[i]]] = PetscMax(PetscAbsScalar(b_aij->a[i]),work[garray[b_aij->j[i]]]); 147 } 148 149 } else SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONG,"Unknown NormType"); 150 if (type == NORM_INFINITY) { 151 ierr = MPI_Allreduce(work,norms,n,MPIU_REAL,MPI_MAX,A->hdr.comm);CHKERRQ(ierr); 152 } else { 153 ierr = MPI_Allreduce(work,norms,n,MPIU_REAL,MPI_SUM,A->hdr.comm);CHKERRQ(ierr); 154 } 155 ierr = PetscFree(work);CHKERRQ(ierr); 156 if (type == NORM_2) { 157 for (i=0; i<n; i++) norms[i] = sqrt(norms[i]); 158 } 159 PetscFunctionReturn(0); 160 } 161 162 #undef __FUNCT__ 163 #define __FUNCT__ "MatGetColumnNorms_SeqDense" 164 PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms) 165 { 166 PetscErrorCode ierr; 167 PetscInt i,j,m,n; 168 PetscScalar *a; 169 170 PetscFunctionBegin; 171 ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr); 172 ierr = PetscMemzero(norms,n*sizeof(PetscReal));CHKERRQ(ierr); 173 ierr = MatGetArray(A,&a);CHKERRQ(ierr); 174 if (type == NORM_2) { 175 for (i=0; i<n; i++ ){ 176 for (j=0; j<m; j++) { 177 norms[i] += PetscAbsScalar(a[j]*a[j]); 178 } 179 a += m; 180 } 181 } else if (type == NORM_1) { 182 for (i=0; i<n; i++ ){ 183 for (j=0; j<m; j++) { 184 norms[i] += PetscAbsScalar(a[j]); 185 } 186 a += m; 187 } 188 } else if (type == NORM_INFINITY) { 189 for (i=0; i<n; i++ ){ 190 for (j=0; j<m; j++) { 191 norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]); 192 } 193 a += m; 194 } 195 } else SETERRQ(((PetscObject)A)->comm,PETSC_ERR_ARG_WRONG,"Unknown NormType"); 196 if (type == NORM_2) { 197 for (i=0; i<n; i++) norms[i] = sqrt(norms[i]); 198 } 199 PetscFunctionReturn(0); 200 } 201 202 #include "../src/mat/impls/dense/mpi/mpidense.h" 203 #undef __FUNCT__ 204 #define __FUNCT__ "MatGetColumnNorms_MPIDense" 205 PetscErrorCode MatGetColumnNorms_MPIDense(Mat A,NormType type,PetscReal *norms) 206 { 207 PetscErrorCode ierr; 208 PetscInt i,n; 209 Mat_MPIDense *a = (Mat_MPIDense*) A->data; 210 PetscReal *work; 211 212 PetscFunctionBegin; 213 ierr = MatGetSize(A,PETSC_NULL,&n);CHKERRQ(ierr); 214 ierr = PetscMalloc(n*sizeof(PetscReal),&work);CHKERRQ(ierr); 215 ierr = MatGetColumnNorms_SeqDense(a->A,type,work);CHKERRQ(ierr); 216 if (type == NORM_2) { 217 for (i=0; i<n; i++) work[i] *= work[i]; 218 } 219 if (type == NORM_INFINITY) { 220 ierr = MPI_Allreduce(work,norms,n,MPIU_REAL,MPI_MAX,A->hdr.comm);CHKERRQ(ierr); 221 } else { 222 ierr = MPI_Allreduce(work,norms,n,MPIU_REAL,MPI_SUM,A->hdr.comm);CHKERRQ(ierr); 223 } 224 ierr = PetscFree(work);CHKERRQ(ierr); 225 if (type == NORM_2) { 226 for (i=0; i<n; i++) norms[i] = sqrt(norms[i]); 227 } 228 PetscFunctionReturn(0); 229 } 230 231 #undef __FUNCT__ 232 #define __FUNCT__ "MatGetColumnNorms" 233 /*@ 234 MatGetColumnNorms - Gets the 2 norms of each column of a sparse or dense matrix. 235 236 Input Parameter: 237 + A - the matrix 238 - type - NORM_2, NORM_1 or NORM_INFINITY 239 240 Output Parameter: 241 . norms - an array as large as the TOTAL number of columns in the matrix 242 243 Level: intermediate 244 245 Notes: Each process has ALL the column norms after the call. Because of the way this is computed each process gets all the values, 246 if each process wants only some of the values it should extract the ones it wants from the array. 247 248 .seealso: MatGetColumns() 249 250 @*/ 251 PetscErrorCode MatGetColumnNorms(Mat A,NormType type,PetscReal *norms) 252 { 253 PetscErrorCode ierr; 254 PetscBool flg; 255 256 PetscFunctionBegin; 257 ierr = PetscTypeCompare((PetscObject)A,MATSEQAIJ,&flg);CHKERRQ(ierr); 258 if (flg) { 259 ierr = MatGetColumnNorms_SeqAIJ(A,type,norms);CHKERRQ(ierr); 260 } else { 261 ierr = PetscTypeCompare((PetscObject)A,MATSEQDENSE,&flg);CHKERRQ(ierr); 262 if (flg) { 263 ierr = MatGetColumnNorms_SeqDense(A,type,norms);CHKERRQ(ierr); 264 } else { 265 ierr = PetscTypeCompare((PetscObject)A,MATMPIDENSE,&flg);CHKERRQ(ierr); 266 if (flg) { 267 ierr = MatGetColumnNorms_MPIDense(A,type,norms);CHKERRQ(ierr); 268 } else { 269 ierr = PetscTypeCompare((PetscObject)A,MATMPIAIJ,&flg);CHKERRQ(ierr); 270 if (flg) { 271 ierr = MatGetColumnNorms_MPIAIJ(A,type,norms);CHKERRQ(ierr); 272 } else SETERRQ(((PetscObject)A)->comm,PETSC_ERR_SUP,"Not coded for this matrix type"); 273 } 274 } 275 } 276 PetscFunctionReturn(0); 277 } 278 279