1 #define PETSCMAT_DLL 2 3 #include "private/matimpl.h" /*I "petscmat.h" I*/ 4 5 #undef __FUNCT__ 6 #define __FUNCT__ "MatGetColumnVector" 7 /*@ 8 MatGetColumnVector - Gets the values from a given column of a matrix. 9 10 Not Collective 11 12 Input Parameters: 13 + A - the matrix 14 . yy - the vector 15 - c - the column requested (in global numbering) 16 17 Level: advanced 18 19 Notes: 20 Each processor for which this is called gets the values for its rows. 21 22 Since PETSc matrices are usually stored in compressed row format, this routine 23 will generally be slow. 24 25 The vector must have the same parallel row layout as the matrix. 26 27 Contributed by: Denis Vanderstraeten 28 29 .keywords: matrix, column, get 30 31 .seealso: MatGetRow(), MatGetDiagonal() 32 33 @*/ 34 PetscErrorCode PETSCMAT_DLLEXPORT MatGetColumnVector(Mat A,Vec yy,PetscInt col) 35 { 36 PetscScalar *y; 37 const PetscScalar *v; 38 PetscErrorCode ierr; 39 PetscInt i,j,nz,N,Rs,Re,rs,re; 40 const PetscInt *idx; 41 42 PetscFunctionBegin; 43 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 44 PetscValidHeaderSpecific(yy,VEC_CLASSID,2); 45 if (col < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Requested negative column: %D",col); 46 ierr = MatGetSize(A,PETSC_NULL,&N);CHKERRQ(ierr); 47 if (col >= N) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Requested column %D larger than number columns in matrix %D",col,N); 48 ierr = MatGetOwnershipRange(A,&Rs,&Re);CHKERRQ(ierr); 49 ierr = VecGetOwnershipRange(yy,&rs,&re);CHKERRQ(ierr); 50 if (Rs != rs || Re != re) SETERRQ4(PETSC_ERR_ARG_INCOMP,"Matrix %D %D does not have same ownership range (size) as vector %D %D",Rs,Re,rs,re); 51 52 if (A->ops->getcolumnvector) { 53 ierr = (*A->ops->getcolumnvector)(A,yy,col);CHKERRQ(ierr); 54 } else { 55 ierr = VecSet(yy,0.0);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 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 74 } 75 PetscFunctionReturn(0); 76 } 77 78 #include "../src/mat/impls/aij/seq/aij.h" 79 80 #undef __FUNCT__ 81 #define __FUNCT__ "MatGetColumnNorms_SeqAIJ" 82 PetscErrorCode MatGetColumnNorms_SeqAIJ(Mat A,NormType type,PetscReal *norms) 83 { 84 PetscErrorCode ierr; 85 PetscInt i,m,n; 86 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 87 88 PetscFunctionBegin; 89 ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr); 90 ierr = PetscMemzero(norms,n*sizeof(PetscReal));CHKERRQ(ierr); 91 if (type == NORM_2) { 92 for (i=0; i<aij->i[m]; i++) { 93 norms[aij->j[i]] += PetscAbsScalar(aij->a[i]*aij->a[i]); 94 } 95 } else if (type == NORM_1) { 96 for (i=0; i<aij->i[m]; i++) { 97 norms[aij->j[i]] += PetscAbsScalar(aij->a[i]); 98 } 99 } else if (type == NORM_INFINITY) { 100 for (i=0; i<aij->i[m]; i++) { 101 norms[aij->j[i]] = PetscMax(PetscAbsScalar(aij->a[i]),norms[aij->j[i]]); 102 } 103 } else SETERRQ(PETSC_ERR_ARG_WRONG,"Unknown NormType"); 104 105 if (type == NORM_2) { 106 for (i=0; i<n; i++) norms[i] = sqrt(norms[i]); 107 } 108 PetscFunctionReturn(0); 109 } 110 111 #include "../src/mat/impls/aij/mpi/mpiaij.h" 112 113 #undef __FUNCT__ 114 #define __FUNCT__ "MatGetColumnNorms_MPIAIJ" 115 PetscErrorCode MatGetColumnNorms_MPIAIJ(Mat A,NormType type,PetscReal *norms) 116 { 117 PetscErrorCode ierr; 118 Mat_MPIAIJ *aij = (Mat_MPIAIJ*)A->data; 119 PetscInt i,n,*garray = aij->garray; 120 Mat_SeqAIJ *a_aij = (Mat_SeqAIJ*) aij->A->data; 121 Mat_SeqAIJ *b_aij = (Mat_SeqAIJ*) aij->B->data; 122 PetscReal *work; 123 124 PetscFunctionBegin; 125 ierr = MatGetSize(A,PETSC_NULL,&n);CHKERRQ(ierr); 126 ierr = PetscMalloc(n*sizeof(PetscReal),&work);CHKERRQ(ierr); 127 ierr = PetscMemzero(work,n*sizeof(PetscReal));CHKERRQ(ierr); 128 if (type == NORM_2) { 129 for (i=0; i<a_aij->i[aij->A->rmap->n]; i++) { 130 work[A->rmap->rstart + a_aij->j[i]] += PetscAbsScalar(a_aij->a[i]*a_aij->a[i]); 131 } 132 for (i=0; i<b_aij->i[aij->B->rmap->n]; i++) { 133 work[garray[b_aij->j[i]]] += PetscAbsScalar(b_aij->a[i]*b_aij->a[i]); 134 } 135 } else if (type == NORM_1) { 136 for (i=0; i<a_aij->i[aij->A->rmap->n]; i++) { 137 work[A->rmap->rstart + a_aij->j[i]] += PetscAbsScalar(a_aij->a[i]); 138 } 139 for (i=0; i<b_aij->i[aij->B->rmap->n]; i++) { 140 work[garray[b_aij->j[i]]] += PetscAbsScalar(b_aij->a[i]); 141 } 142 } else if (type == NORM_INFINITY) { 143 for (i=0; i<a_aij->i[aij->A->rmap->n]; i++) { 144 work[A->rmap->rstart + a_aij->j[i]] = PetscMax(PetscAbsScalar(a_aij->a[i]), work[A->rmap->rstart + a_aij->j[i]]); 145 } 146 for (i=0; i<b_aij->i[aij->B->rmap->n]; i++) { 147 work[garray[b_aij->j[i]]] = PetscMax(PetscAbsScalar(b_aij->a[i]),work[garray[b_aij->j[i]]]); 148 } 149 150 } else SETERRQ(PETSC_ERR_ARG_WRONG,"Unknown NormType"); 151 if (type == NORM_INFINITY) { 152 ierr = MPI_Allreduce(work,norms,n,MPIU_REAL,MPI_MAX,A->hdr.comm);CHKERRQ(ierr); 153 } else { 154 ierr = MPI_Allreduce(work,norms,n,MPIU_REAL,MPI_SUM,A->hdr.comm);CHKERRQ(ierr); 155 } 156 ierr = PetscFree(work);CHKERRQ(ierr); 157 if (type == NORM_2) { 158 for (i=0; i<n; i++) norms[i] = sqrt(norms[i]); 159 } 160 PetscFunctionReturn(0); 161 } 162 163 #undef __FUNCT__ 164 #define __FUNCT__ "MatGetColumnNorms_SeqDense" 165 PetscErrorCode MatGetColumnNorms_SeqDense(Mat A,NormType type,PetscReal *norms) 166 { 167 PetscErrorCode ierr; 168 PetscInt i,j,m,n; 169 PetscScalar *a; 170 171 PetscFunctionBegin; 172 ierr = MatGetSize(A,&m,&n);CHKERRQ(ierr); 173 ierr = PetscMemzero(norms,n*sizeof(PetscReal));CHKERRQ(ierr); 174 ierr = MatGetArray(A,&a);CHKERRQ(ierr); 175 if (type == NORM_2) { 176 for (i=0; i<n; i++ ){ 177 for (j=0; j<m; j++) { 178 norms[i] += PetscAbsScalar(a[j]*a[j]); 179 } 180 a += m; 181 } 182 } else if (type == NORM_1) { 183 for (i=0; i<n; i++ ){ 184 for (j=0; j<m; j++) { 185 norms[i] += PetscAbsScalar(a[j]); 186 } 187 a += m; 188 } 189 } else if (type == NORM_INFINITY) { 190 for (i=0; i<n; i++ ){ 191 for (j=0; j<m; j++) { 192 norms[i] = PetscMax(PetscAbsScalar(a[j]),norms[i]); 193 } 194 a += m; 195 } 196 } else SETERRQ(PETSC_ERR_ARG_WRONG,"Unknown NormType"); 197 if (type == NORM_2) { 198 for (i=0; i<n; i++) norms[i] = sqrt(norms[i]); 199 } 200 PetscFunctionReturn(0); 201 } 202 203 #include "../src/mat/impls/dense/mpi/mpidense.h" 204 #undef __FUNCT__ 205 #define __FUNCT__ "MatGetColumnNorms_MPIDense" 206 PetscErrorCode MatGetColumnNorms_MPIDense(Mat A,NormType type,PetscReal *norms) 207 { 208 PetscErrorCode ierr; 209 PetscInt i,n; 210 Mat_MPIDense *a = (Mat_MPIDense*) A->data; 211 PetscReal *work; 212 213 PetscFunctionBegin; 214 ierr = MatGetSize(A,PETSC_NULL,&n);CHKERRQ(ierr); 215 ierr = PetscMalloc(n*sizeof(PetscReal),&work);CHKERRQ(ierr); 216 ierr = MatGetColumnNorms_SeqDense(a->A,type,work);CHKERRQ(ierr); 217 if (type == NORM_2) { 218 for (i=0; i<n; i++) work[i] *= work[i]; 219 } 220 if (type == NORM_INFINITY) { 221 ierr = MPI_Allreduce(work,norms,n,MPIU_REAL,MPI_MAX,A->hdr.comm);CHKERRQ(ierr); 222 } else { 223 ierr = MPI_Allreduce(work,norms,n,MPIU_REAL,MPI_SUM,A->hdr.comm);CHKERRQ(ierr); 224 } 225 ierr = PetscFree(work);CHKERRQ(ierr); 226 if (type == NORM_2) { 227 for (i=0; i<n; i++) norms[i] = sqrt(norms[i]); 228 } 229 PetscFunctionReturn(0); 230 } 231 232 #undef __FUNCT__ 233 #define __FUNCT__ "MatGetColumnNorms" 234 /*@ 235 MatGetColumnNorms - Gets the 2 norms of each column of a sparse or dense matrix. 236 237 Input Parameter: 238 + A - the matrix 239 - type - NORM_2, NORM_1 or NORM_INFINITY 240 241 Output Parameter: 242 . norms - an array as large as the TOTAL number of columns in the matrix 243 244 Notes: Each process has ALL the column norms after the call. Because of the way this is computed each process gets all the values, 245 if each process wants only some of the values it should extract the ones it wants from the array. 246 247 .seealso: MatGetColumns() 248 249 @*/ 250 PetscErrorCode MatGetColumnNorms(Mat A,NormType type,PetscReal *norms) 251 { 252 PetscErrorCode ierr; 253 PetscTruth flg; 254 255 PetscFunctionBegin; 256 ierr = PetscTypeCompare((PetscObject)A,MATSEQAIJ,&flg);CHKERRQ(ierr); 257 if (flg) { 258 ierr = MatGetColumnNorms_SeqAIJ(A,type,norms);CHKERRQ(ierr); 259 } else { 260 ierr = PetscTypeCompare((PetscObject)A,MATSEQDENSE,&flg);CHKERRQ(ierr); 261 if (flg) { 262 ierr = MatGetColumnNorms_SeqDense(A,type,norms);CHKERRQ(ierr); 263 } else { 264 ierr = PetscTypeCompare((PetscObject)A,MATMPIDENSE,&flg);CHKERRQ(ierr); 265 if (flg) { 266 ierr = MatGetColumnNorms_MPIDense(A,type,norms);CHKERRQ(ierr); 267 } else { 268 ierr = PetscTypeCompare((PetscObject)A,MATMPIAIJ,&flg);CHKERRQ(ierr); 269 if (flg) { 270 ierr = MatGetColumnNorms_MPIAIJ(A,type,norms);CHKERRQ(ierr); 271 } else SETERRQ(PETSC_ERR_SUP,"Not coded for this matrix type"); 272 } 273 } 274 } 275 PetscFunctionReturn(0); 276 } 277 278