1 2 #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/ 3 4 /*@ 5 MatGetColumnVector - Gets the values from a given column of a matrix. 6 7 Not Collective 8 9 Input Parameters: 10 + A - the matrix 11 . yy - the vector 12 - col - the column requested (in global numbering) 13 14 Level: advanced 15 16 Notes: 17 If a Mat type does not implement the operation, each processor for which this is called 18 gets the values for its rows using MatGetRow(). 19 20 The vector must have the same parallel row layout as the matrix. 21 22 Contributed by: Denis Vanderstraeten 23 24 .seealso: MatGetRow(), MatGetDiagonal(), MatMult() 25 26 @*/ 27 PetscErrorCode MatGetColumnVector(Mat A,Vec yy,PetscInt col) 28 { 29 PetscScalar *y; 30 const PetscScalar *v; 31 PetscErrorCode ierr; 32 PetscInt i,j,nz,N,Rs,Re,rs,re; 33 const PetscInt *idx; 34 35 PetscFunctionBegin; 36 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 37 PetscValidHeaderSpecific(yy,VEC_CLASSID,2); 38 PetscValidLogicalCollectiveInt(A,col,3); 39 if (col < 0) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Requested negative column: %D",col); 40 ierr = MatGetSize(A,NULL,&N);CHKERRQ(ierr); 41 if (col >= N) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Requested column %D larger than number columns in matrix %D",col,N); 42 ierr = MatGetOwnershipRange(A,&Rs,&Re);CHKERRQ(ierr); 43 ierr = VecGetOwnershipRange(yy,&rs,&re);CHKERRQ(ierr); 44 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); 45 46 if (A->ops->getcolumnvector) { 47 ierr = (*A->ops->getcolumnvector)(A,yy,col);CHKERRQ(ierr); 48 } else { 49 ierr = VecSet(yy,0.0);CHKERRQ(ierr); 50 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 51 /* TODO for general matrices */ 52 for (i=Rs; i<Re; i++) { 53 ierr = MatGetRow(A,i,&nz,&idx,&v);CHKERRQ(ierr); 54 if (nz && idx[0] <= col) { 55 /* 56 Should use faster search here 57 */ 58 for (j=0; j<nz; j++) { 59 if (idx[j] >= col) { 60 if (idx[j] == col) y[i-rs] = v[j]; 61 break; 62 } 63 } 64 } 65 ierr = MatRestoreRow(A,i,&nz,&idx,&v);CHKERRQ(ierr); 66 } 67 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 68 } 69 PetscFunctionReturn(0); 70 } 71 72 /*@ 73 MatGetColumnNorms - Gets the norms of each column of a sparse or dense matrix. 74 75 Input Parameters: 76 + A - the matrix 77 - type - NORM_2, NORM_1 or NORM_INFINITY 78 79 Output Parameter: 80 . norms - an array as large as the TOTAL number of columns in the matrix 81 82 Level: intermediate 83 84 Notes: 85 Each process has ALL the column norms after the call. Because of the way this is computed each process gets all the values, 86 if each process wants only some of the values it should extract the ones it wants from the array. 87 88 .seealso: NormType, MatNorm() 89 90 @*/ 91 PetscErrorCode MatGetColumnNorms(Mat A,NormType type,PetscReal norms[]) 92 { 93 /* NOTE: MatGetColumnNorms() could simply be a macro that calls MatGetColumnReductions(). 94 * I've kept this as a function because it allows slightly more in the way of error checking, 95 * erroring out if MatGetColumnNorms() is not called with a valid NormType. */ 96 PetscErrorCode ierr; 97 98 PetscFunctionBegin; 99 if (type == NORM_2 || type == NORM_1 || type == NORM_FROBENIUS || type == NORM_INFINITY || type == NORM_1_AND_2) { 100 ierr = MatGetColumnReductions(A,type,norms);CHKERRQ(ierr); 101 } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Unknown NormType"); 102 PetscFunctionReturn(0); 103 } 104 105 /*@ 106 MatGetColumnSumsRealPart - Gets the sums of the real part of each column of a sparse or dense matrix. 107 108 Input Parameter: 109 . A - the matrix 110 111 Output Parameter: 112 . sums - an array as large as the TOTAL number of columns in the matrix 113 114 Level: intermediate 115 116 Notes: 117 Each process has ALL the column sums after the call. Because of the way this is computed each process gets all the values, 118 if each process wants only some of the values it should extract the ones it wants from the array. 119 120 .seealso: MatGetColumnSumsImaginaryPart(), VecSum(), MatGetColumnMeans(), MatGetColumnNorms(), MatGetColumnReductions() 121 122 @*/ 123 PetscErrorCode MatGetColumnSumsRealPart(Mat A,PetscReal sums[]) 124 { 125 PetscErrorCode ierr; 126 127 PetscFunctionBegin; 128 ierr = MatGetColumnReductions(A,REDUCTION_SUM_REALPART,sums);CHKERRQ(ierr); 129 PetscFunctionReturn(0); 130 } 131 132 /*@ 133 MatGetColumnSumsImaginaryPart - Gets the sums of the imaginary part of each column of a sparse or dense matrix. 134 135 Input Parameter: 136 . A - the matrix 137 138 Output Parameter: 139 . sums - an array as large as the TOTAL number of columns in the matrix 140 141 Level: intermediate 142 143 Notes: 144 Each process has ALL the column sums after the call. Because of the way this is computed each process gets all the values, 145 if each process wants only some of the values it should extract the ones it wants from the array. 146 147 .seealso: MatGetColumnSumsRealPart(), VecSum(), MatGetColumnMeans(), MatGetColumnNorms(), MatGetColumnReductions() 148 149 @*/ 150 PetscErrorCode MatGetColumnSumsImaginaryPart(Mat A,PetscReal sums[]) 151 { 152 PetscErrorCode ierr; 153 154 PetscFunctionBegin; 155 ierr = MatGetColumnReductions(A,REDUCTION_SUM_IMAGINARYPART,sums);CHKERRQ(ierr); 156 PetscFunctionReturn(0); 157 } 158 159 /*@ 160 MatGetColumnSums - Gets the sums of each column of a sparse or dense matrix. 161 162 Input Parameter: 163 . A - the matrix 164 165 Output Parameter: 166 . sums - an array as large as the TOTAL number of columns in the matrix 167 168 Level: intermediate 169 170 Notes: 171 Each process has ALL the column sums after the call. Because of the way this is computed each process gets all the values, 172 if each process wants only some of the values it should extract the ones it wants from the array. 173 174 .seealso: VecSum(), MatGetColumnMeans(), MatGetColumnNorms(), MatGetColumnReductions() 175 176 @*/ 177 PetscErrorCode MatGetColumnSums(Mat A,PetscScalar sums[]) 178 { 179 PetscErrorCode ierr; 180 #if defined(PETSC_USE_COMPLEX) 181 PetscInt i,n; 182 PetscReal *work; 183 #endif 184 185 PetscFunctionBegin; 186 187 #if !defined(PETSC_USE_COMPLEX) 188 ierr = MatGetColumnSumsRealPart(A,sums);CHKERRQ(ierr); 189 #else 190 ierr = MatGetSize(A,NULL,&n);CHKERRQ(ierr); 191 ierr = PetscArrayzero(sums,n);CHKERRQ(ierr); 192 ierr = PetscCalloc1(n,&work);CHKERRQ(ierr); 193 ierr = MatGetColumnSumsRealPart(A,work);CHKERRQ(ierr); 194 for (i=0; i<n; i++) sums[i] = work[i]; 195 ierr = MatGetColumnSumsImaginaryPart(A,work);CHKERRQ(ierr); 196 for (i=0; i<n; i++) sums[i] += work[i]*PETSC_i; 197 ierr = PetscFree(work);CHKERRQ(ierr); 198 #endif 199 PetscFunctionReturn(0); 200 } 201 202 /*@ 203 MatGetColumnMeansRealPart - Gets the arithmetic means of the real part of each column of a sparse or dense matrix. 204 205 Input Parameter: 206 . A - the matrix 207 208 Output Parameter: 209 . sums - an array as large as the TOTAL number of columns in the matrix 210 211 Level: intermediate 212 213 Notes: 214 Each process has ALL the column means after the call. Because of the way this is computed each process gets all the values, 215 if each process wants only some of the values it should extract the ones it wants from the array. 216 217 .seealso: MatGetColumnMeansImaginaryPart(), VecSum(), MatGetColumnSums(), MatGetColumnNorms(), MatGetColumnReductions() 218 219 @*/ 220 PetscErrorCode MatGetColumnMeansRealPart(Mat A,PetscReal means[]) 221 { 222 PetscErrorCode ierr; 223 224 PetscFunctionBegin; 225 ierr = MatGetColumnReductions(A,REDUCTION_MEAN_REALPART,means);CHKERRQ(ierr); 226 PetscFunctionReturn(0); 227 } 228 229 /*@ 230 MatGetColumnMeansImaginaryPart - Gets the arithmetic means of the imaginary part of each column of a sparse or dense matrix. 231 232 Input Parameter: 233 . A - the matrix 234 235 Output Parameter: 236 . sums - an array as large as the TOTAL number of columns in the matrix 237 238 Level: intermediate 239 240 Notes: 241 Each process has ALL the column means after the call. Because of the way this is computed each process gets all the values, 242 if each process wants only some of the values it should extract the ones it wants from the array. 243 244 .seealso: MatGetColumnMeansRealPart(), VecSum(), MatGetColumnSums(), MatGetColumnNorms(), MatGetColumnReductions() 245 246 @*/ 247 PetscErrorCode MatGetColumnMeansImaginaryPart(Mat A,PetscReal means[]) 248 { 249 PetscErrorCode ierr; 250 251 PetscFunctionBegin; 252 ierr = MatGetColumnReductions(A,REDUCTION_MEAN_IMAGINARYPART,means);CHKERRQ(ierr); 253 PetscFunctionReturn(0); 254 } 255 256 /*@ 257 MatGetColumnMeans - Gets the arithmetic means of each column of a sparse or dense matrix. 258 259 Input Parameter: 260 . A - the matrix 261 262 Output Parameter: 263 . means - an array as large as the TOTAL number of columns in the matrix 264 265 Level: intermediate 266 267 Notes: 268 Each process has ALL the column means after the call. Because of the way this is computed each process gets all the values, 269 if each process wants only some of the values it should extract the ones it wants from the array. 270 271 .seealso: VecSum(), MatGetColumnSums(), MatGetColumnNorms(), MatGetColumnReductions() 272 273 @*/ 274 PetscErrorCode MatGetColumnMeans(Mat A,PetscScalar means[]) 275 { 276 PetscErrorCode ierr; 277 #if defined(PETSC_USE_COMPLEX) 278 PetscInt i,n; 279 PetscReal *work; 280 #endif 281 282 PetscFunctionBegin; 283 284 #if !defined(PETSC_USE_COMPLEX) 285 ierr = MatGetColumnMeansRealPart(A,means);CHKERRQ(ierr); 286 #else 287 ierr = MatGetSize(A,NULL,&n);CHKERRQ(ierr); 288 ierr = PetscArrayzero(means,n);CHKERRQ(ierr); 289 ierr = PetscCalloc1(n,&work);CHKERRQ(ierr); 290 ierr = MatGetColumnMeansRealPart(A,work);CHKERRQ(ierr); 291 for (i=0; i<n; i++) means[i] = work[i]; 292 ierr = MatGetColumnMeansImaginaryPart(A,work);CHKERRQ(ierr); 293 for (i=0; i<n; i++) means[i] += work[i]*PETSC_i; 294 ierr = PetscFree(work);CHKERRQ(ierr); 295 #endif 296 PetscFunctionReturn(0); 297 } 298 299 /*@ 300 MatGetColumnReductions - Gets the reductions of each column of a sparse or dense matrix. 301 302 Input Parameters: 303 + A - the matrix 304 - type - A constant defined in NormType or ReductionType: NORM_2, NORM_1, NORM_INFINITY, REDUCTION_SUM_REALPART, 305 REDUCTION_SUM_IMAGINARYPART, REDUCTION_MEAN_REALPART, REDUCTION_MEAN_IMAGINARYPART 306 307 Output Parameter: 308 . reductions - an array as large as the TOTAL number of columns in the matrix 309 310 Level: developer 311 312 Notes: 313 Each process has ALL the column reductions after the call. Because of the way this is computed each process gets all the values, 314 if each process wants only some of the values it should extract the ones it wants from the array. 315 316 Developer Note: 317 This routine is primarily intended as a back-end. 318 MatGetColumnNorms(), MatGetColumnSums(), and MatGetColumnMeans() are implemented using this routine. 319 320 .seealso: ReductionType, NormType, MatGetColumnNorms(), MatGetColumnSums(), MatGetColumnMeans() 321 322 @*/ 323 PetscErrorCode MatGetColumnReductions(Mat A,PetscInt type,PetscReal reductions[]) 324 { 325 PetscErrorCode ierr; 326 327 PetscFunctionBegin; 328 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 329 if (A->ops->getcolumnreductions) { 330 ierr = (*A->ops->getcolumnreductions)(A,type,reductions);CHKERRQ(ierr); 331 } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Not coded for this matrix type"); 332 PetscFunctionReturn(0); 333 } 334