1 #define PETSCMAT_DLL 2 3 /* 4 Defines the basic matrix operations for sequential dense. 5 */ 6 7 #include "src/mat/impls/dense/seq/dense.h" 8 #include "petscblaslapack.h" 9 10 #undef __FUNCT__ 11 #define __FUNCT__ "MatAXPY_SeqDense" 12 PetscErrorCode MatAXPY_SeqDense(Mat Y,PetscScalar alpha,Mat X,MatStructure str) 13 { 14 Mat_SeqDense *x = (Mat_SeqDense*)X->data,*y = (Mat_SeqDense*)Y->data; 15 PetscScalar oalpha = alpha; 16 PetscInt j; 17 PetscBLASInt N,m,ldax,lday,one = 1; 18 PetscErrorCode ierr; 19 20 PetscFunctionBegin; 21 N = PetscBLASIntCast(X->rmap->n*X->cmap->n); 22 m = PetscBLASIntCast(X->rmap->n); 23 ldax = PetscBLASIntCast(x->lda); 24 lday = PetscBLASIntCast(y->lda); 25 if (ldax>m || lday>m) { 26 for (j=0; j<X->cmap->n; j++) { 27 BLASaxpy_(&m,&oalpha,x->v+j*ldax,&one,y->v+j*lday,&one); 28 } 29 } else { 30 BLASaxpy_(&N,&oalpha,x->v,&one,y->v,&one); 31 } 32 ierr = PetscLogFlops(2*N-1);CHKERRQ(ierr); 33 PetscFunctionReturn(0); 34 } 35 36 #undef __FUNCT__ 37 #define __FUNCT__ "MatGetInfo_SeqDense" 38 PetscErrorCode MatGetInfo_SeqDense(Mat A,MatInfoType flag,MatInfo *info) 39 { 40 PetscInt N = A->rmap->n*A->cmap->n; 41 42 PetscFunctionBegin; 43 info->rows_global = (double)A->rmap->n; 44 info->columns_global = (double)A->cmap->n; 45 info->rows_local = (double)A->rmap->n; 46 info->columns_local = (double)A->cmap->n; 47 info->block_size = 1.0; 48 info->nz_allocated = (double)N; 49 info->nz_used = (double)N; 50 info->nz_unneeded = (double)0; 51 info->assemblies = (double)A->num_ass; 52 info->mallocs = 0; 53 info->memory = ((PetscObject)A)->mem; 54 info->fill_ratio_given = 0; 55 info->fill_ratio_needed = 0; 56 info->factor_mallocs = 0; 57 PetscFunctionReturn(0); 58 } 59 60 #undef __FUNCT__ 61 #define __FUNCT__ "MatScale_SeqDense" 62 PetscErrorCode MatScale_SeqDense(Mat A,PetscScalar alpha) 63 { 64 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 65 PetscScalar oalpha = alpha; 66 PetscErrorCode ierr; 67 PetscBLASInt one = 1,j,nz,lda = PetscBLASIntCast(a->lda); 68 69 PetscFunctionBegin; 70 if (lda>A->rmap->n) { 71 nz = PetscBLASIntCast(A->rmap->n); 72 for (j=0; j<A->cmap->n; j++) { 73 BLASscal_(&nz,&oalpha,a->v+j*lda,&one); 74 } 75 } else { 76 nz = PetscBLASIntCast(A->rmap->n*A->cmap->n); 77 BLASscal_(&nz,&oalpha,a->v,&one); 78 } 79 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 80 PetscFunctionReturn(0); 81 } 82 83 #undef __FUNCT__ 84 #define __FUNCT__ "MatIsHermitian_SeqDense" 85 PetscErrorCode MatIsHermitian_SeqDense(Mat A,PetscReal rtol,PetscTruth *fl) 86 { 87 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 88 PetscInt i,j,m = A->rmap->n,N; 89 PetscScalar *v = a->v; 90 91 PetscFunctionBegin; 92 *fl = PETSC_FALSE; 93 if (A->rmap->n != A->cmap->n) PetscFunctionReturn(0); 94 N = a->lda; 95 96 for (i=0; i<m; i++) { 97 for (j=i+1; j<m; j++) { 98 if (PetscAbsScalar(v[i+j*N] - PetscConj(v[j+i*N])) > rtol) PetscFunctionReturn(0); 99 } 100 } 101 *fl = PETSC_TRUE; 102 PetscFunctionReturn(0); 103 } 104 105 #undef __FUNCT__ 106 #define __FUNCT__ "MatDuplicateNoCreate_SeqDense" 107 PetscErrorCode MatDuplicateNoCreate_SeqDense(Mat A,MatDuplicateOption cpvalues,Mat *newmat) 108 { 109 Mat_SeqDense *mat = (Mat_SeqDense*)A->data,*l; 110 PetscErrorCode ierr; 111 PetscInt lda = (PetscInt)mat->lda,j,m; 112 Mat newi = *newmat; 113 114 PetscFunctionBegin; 115 ierr = MatSeqDenseSetPreallocation(*newmat,PETSC_NULL);CHKERRQ(ierr); 116 if (cpvalues == MAT_COPY_VALUES) { 117 l = (Mat_SeqDense*)newi->data; 118 if (lda>A->rmap->n) { 119 m = A->rmap->n; 120 for (j=0; j<A->cmap->n; j++) { 121 ierr = PetscMemcpy(l->v+j*m,mat->v+j*lda,m*sizeof(PetscScalar));CHKERRQ(ierr); 122 } 123 } else { 124 ierr = PetscMemcpy(l->v,mat->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 125 } 126 } 127 newi->assembled = PETSC_TRUE; 128 PetscFunctionReturn(0); 129 } 130 131 #undef __FUNCT__ 132 #define __FUNCT__ "MatDuplicate_SeqDense" 133 PetscErrorCode MatDuplicate_SeqDense(Mat A,MatDuplicateOption cpvalues,Mat *newmat) 134 { 135 PetscErrorCode ierr; 136 137 PetscFunctionBegin; 138 ierr = MatCreate(((PetscObject)A)->comm,newmat);CHKERRQ(ierr); 139 ierr = MatSetSizes(*newmat,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 140 ierr = MatSetType(*newmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 141 ierr = MatDuplicateNoCreate_SeqDense(A,cpvalues,newmat);CHKERRQ(ierr); 142 PetscFunctionReturn(0); 143 } 144 145 EXTERN PetscErrorCode MatLUFactor_SeqDense(Mat A,IS row,IS col,MatFactorInfo *minfo); 146 147 #undef __FUNCT__ 148 #define __FUNCT__ "MatLUFactorNumeric_SeqDense" 149 PetscErrorCode MatLUFactorNumeric_SeqDense(Mat A,MatFactorInfo *info_dummy,Mat *fact) 150 { 151 PetscErrorCode ierr; 152 MatFactorInfo info; 153 154 PetscFunctionBegin; 155 /* copy the numerical values */ 156 ierr = MatDuplicateNoCreate_SeqDense(A,MAT_COPY_VALUES,fact);CHKERRQ(ierr); 157 ierr = MatLUFactor_SeqDense(*fact,0,0,&info);CHKERRQ(ierr) 158 PetscFunctionReturn(0); 159 } 160 161 #undef __FUNCT__ 162 #define __FUNCT__ "MatSolve_SeqDense" 163 PetscErrorCode MatSolve_SeqDense(Mat A,Vec xx,Vec yy) 164 { 165 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 166 PetscErrorCode ierr; 167 PetscScalar *x,*y; 168 PetscBLASInt one = 1,info,m = PetscBLASIntCast(A->rmap->n); 169 170 PetscFunctionBegin; 171 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 172 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 173 ierr = PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 174 if (A->factor == MAT_FACTOR_LU) { 175 #if defined(PETSC_MISSING_LAPACK_GETRS) 176 SETERRQ(PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable."); 177 #else 178 LAPACKgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info); 179 if (info) SETERRQ(PETSC_ERR_LIB,"GETRS - Bad solve"); 180 #endif 181 } else if (A->factor == MAT_FACTOR_CHOLESKY){ 182 #if defined(PETSC_MISSING_LAPACK_POTRS) 183 SETERRQ(PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable."); 184 #else 185 LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info); 186 if (info) SETERRQ(PETSC_ERR_LIB,"POTRS Bad solve"); 187 #endif 188 } 189 else SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Matrix must be factored to solve"); 190 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 191 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 192 ierr = PetscLogFlops(2*A->cmap->n*A->cmap->n - A->cmap->n);CHKERRQ(ierr); 193 PetscFunctionReturn(0); 194 } 195 196 #undef __FUNCT__ 197 #define __FUNCT__ "MatSolveTranspose_SeqDense" 198 PetscErrorCode MatSolveTranspose_SeqDense(Mat A,Vec xx,Vec yy) 199 { 200 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 201 PetscErrorCode ierr; 202 PetscScalar *x,*y; 203 PetscBLASInt one = 1,info,m = PetscBLASIntCast(A->rmap->n); 204 205 PetscFunctionBegin; 206 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 207 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 208 ierr = PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 209 /* assume if pivots exist then use LU; else Cholesky */ 210 if (mat->pivots) { 211 #if defined(PETSC_MISSING_LAPACK_GETRS) 212 SETERRQ(PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable."); 213 #else 214 LAPACKgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info); 215 if (info) SETERRQ(PETSC_ERR_LIB,"POTRS - Bad solve"); 216 #endif 217 } else { 218 #if defined(PETSC_MISSING_LAPACK_POTRS) 219 SETERRQ(PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable."); 220 #else 221 LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info); 222 if (info) SETERRQ(PETSC_ERR_LIB,"POTRS - Bad solve"); 223 #endif 224 } 225 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 226 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 227 ierr = PetscLogFlops(2*A->cmap->n*A->cmap->n - A->cmap->n);CHKERRQ(ierr); 228 PetscFunctionReturn(0); 229 } 230 231 #undef __FUNCT__ 232 #define __FUNCT__ "MatSolveAdd_SeqDense" 233 PetscErrorCode MatSolveAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 234 { 235 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 236 PetscErrorCode ierr; 237 PetscScalar *x,*y,sone = 1.0; 238 Vec tmp = 0; 239 PetscBLASInt one = 1,info,m = PetscBLASIntCast(A->rmap->n); 240 241 PetscFunctionBegin; 242 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 243 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 244 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 245 if (yy == zz) { 246 ierr = VecDuplicate(yy,&tmp);CHKERRQ(ierr); 247 ierr = PetscLogObjectParent(A,tmp);CHKERRQ(ierr); 248 ierr = VecCopy(yy,tmp);CHKERRQ(ierr); 249 } 250 ierr = PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 251 /* assume if pivots exist then use LU; else Cholesky */ 252 if (mat->pivots) { 253 #if defined(PETSC_MISSING_LAPACK_GETRS) 254 SETERRQ(PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable."); 255 #else 256 LAPACKgetrs_("N",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info); 257 if (info) SETERRQ(PETSC_ERR_LIB,"Bad solve"); 258 #endif 259 } else { 260 #if defined(PETSC_MISSING_LAPACK_POTRS) 261 SETERRQ(PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable."); 262 #else 263 LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info); 264 if (info) SETERRQ(PETSC_ERR_LIB,"Bad solve"); 265 #endif 266 } 267 if (tmp) {ierr = VecAXPY(yy,sone,tmp);CHKERRQ(ierr); ierr = VecDestroy(tmp);CHKERRQ(ierr);} 268 else {ierr = VecAXPY(yy,sone,zz);CHKERRQ(ierr);} 269 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 270 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 271 ierr = PetscLogFlops(2*A->cmap->n*A->cmap->n);CHKERRQ(ierr); 272 PetscFunctionReturn(0); 273 } 274 275 #undef __FUNCT__ 276 #define __FUNCT__ "MatSolveTransposeAdd_SeqDense" 277 PetscErrorCode MatSolveTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 278 { 279 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 280 PetscErrorCode ierr; 281 PetscScalar *x,*y,sone = 1.0; 282 Vec tmp; 283 PetscBLASInt one = 1,info,m = PetscBLASIntCast(A->rmap->n); 284 285 PetscFunctionBegin; 286 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 287 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 288 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 289 if (yy == zz) { 290 ierr = VecDuplicate(yy,&tmp);CHKERRQ(ierr); 291 ierr = PetscLogObjectParent(A,tmp);CHKERRQ(ierr); 292 ierr = VecCopy(yy,tmp);CHKERRQ(ierr); 293 } 294 ierr = PetscMemcpy(y,x,A->rmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 295 /* assume if pivots exist then use LU; else Cholesky */ 296 if (mat->pivots) { 297 #if defined(PETSC_MISSING_LAPACK_GETRS) 298 SETERRQ(PETSC_ERR_SUP,"GETRS - Lapack routine is unavailable."); 299 #else 300 LAPACKgetrs_("T",&m,&one,mat->v,&mat->lda,mat->pivots,y,&m,&info); 301 if (info) SETERRQ(PETSC_ERR_LIB,"Bad solve"); 302 #endif 303 } else { 304 #if defined(PETSC_MISSING_LAPACK_POTRS) 305 SETERRQ(PETSC_ERR_SUP,"POTRS - Lapack routine is unavailable."); 306 #else 307 LAPACKpotrs_("L",&m,&one,mat->v,&mat->lda,y,&m,&info); 308 if (info) SETERRQ(PETSC_ERR_LIB,"Bad solve"); 309 #endif 310 } 311 if (tmp) { 312 ierr = VecAXPY(yy,sone,tmp);CHKERRQ(ierr); 313 ierr = VecDestroy(tmp);CHKERRQ(ierr); 314 } else { 315 ierr = VecAXPY(yy,sone,zz);CHKERRQ(ierr); 316 } 317 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 318 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 319 ierr = PetscLogFlops(2*A->cmap->n*A->cmap->n);CHKERRQ(ierr); 320 PetscFunctionReturn(0); 321 } 322 323 /* ---------------------------------------------------------------*/ 324 /* COMMENT: I have chosen to hide row permutation in the pivots, 325 rather than put it in the Mat->row slot.*/ 326 #undef __FUNCT__ 327 #define __FUNCT__ "MatLUFactor_SeqDense" 328 PetscErrorCode MatLUFactor_SeqDense(Mat A,IS row,IS col,MatFactorInfo *minfo) 329 { 330 #if defined(PETSC_MISSING_LAPACK_GETRF) 331 PetscFunctionBegin; 332 SETERRQ(PETSC_ERR_SUP,"GETRF - Lapack routine is unavailable."); 333 #else 334 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 335 PetscErrorCode ierr; 336 PetscBLASInt n,m,info; 337 338 PetscFunctionBegin; 339 n = PetscBLASIntCast(A->cmap->n); 340 m = PetscBLASIntCast(A->rmap->n); 341 if (!mat->pivots) { 342 ierr = PetscMalloc((A->rmap->n+1)*sizeof(PetscBLASInt),&mat->pivots);CHKERRQ(ierr); 343 ierr = PetscLogObjectMemory(A,A->rmap->n*sizeof(PetscBLASInt));CHKERRQ(ierr); 344 } 345 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 346 LAPACKgetrf_(&m,&n,mat->v,&mat->lda,mat->pivots,&info); 347 if (info<0) SETERRQ(PETSC_ERR_LIB,"Bad argument to LU factorization"); 348 if (info>0) SETERRQ(PETSC_ERR_MAT_LU_ZRPVT,"Bad LU factorization"); 349 A->ops->solve = MatSolve_SeqDense; 350 A->ops->solvetranspose = MatSolveTranspose_SeqDense; 351 A->ops->solveadd = MatSolveAdd_SeqDense; 352 A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense; 353 A->factor = MAT_FACTOR_LU; 354 355 ierr = PetscLogFlops((2*A->cmap->n*A->cmap->n*A->cmap->n)/3);CHKERRQ(ierr); 356 #endif 357 PetscFunctionReturn(0); 358 } 359 360 #undef __FUNCT__ 361 #define __FUNCT__ "MatCholeskyFactor_SeqDense" 362 PetscErrorCode MatCholeskyFactor_SeqDense(Mat A,IS perm,MatFactorInfo *factinfo) 363 { 364 #if defined(PETSC_MISSING_LAPACK_POTRF) 365 PetscFunctionBegin; 366 SETERRQ(PETSC_ERR_SUP,"POTRF - Lapack routine is unavailable."); 367 #else 368 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 369 PetscErrorCode ierr; 370 PetscBLASInt info,n = PetscBLASIntCast(A->cmap->n); 371 372 PetscFunctionBegin; 373 ierr = PetscFree(mat->pivots);CHKERRQ(ierr); 374 mat->pivots = 0; 375 376 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 377 LAPACKpotrf_("L",&n,mat->v,&mat->lda,&info); 378 if (info) SETERRQ1(PETSC_ERR_MAT_CH_ZRPVT,"Bad factorization: zero pivot in row %D",(PetscInt)info-1); 379 A->ops->solve = MatSolve_SeqDense; 380 A->ops->solvetranspose = MatSolveTranspose_SeqDense; 381 A->ops->solveadd = MatSolveAdd_SeqDense; 382 A->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense; 383 A->factor = MAT_FACTOR_CHOLESKY; 384 ierr = PetscLogFlops((A->cmap->n*A->cmap->n*A->cmap->n)/3);CHKERRQ(ierr); 385 #endif 386 PetscFunctionReturn(0); 387 } 388 389 390 #undef __FUNCT__ 391 #define __FUNCT__ "MatCholeskyFactorNumeric_SeqDense" 392 PetscErrorCode MatCholeskyFactorNumeric_SeqDense(Mat A,MatFactorInfo *info_dummy,Mat *fact) 393 { 394 PetscErrorCode ierr; 395 MatFactorInfo info; 396 397 PetscFunctionBegin; 398 info.fill = 1.0; 399 ierr = MatDuplicateNoCreate_SeqDense(A,MAT_COPY_VALUES,fact);CHKERRQ(ierr); 400 ierr = MatCholeskyFactor_SeqDense(*fact,0,&info);CHKERRQ(ierr); 401 PetscFunctionReturn(0); 402 } 403 404 #undef __FUNCT__ 405 #define __FUNCT__ "MatCholeskyFactorSymbolic_SeqDense" 406 PetscErrorCode MatCholeskyFactorSymbolic_SeqDense(Mat A,IS row,MatFactorInfo *info,Mat *fact) 407 { 408 PetscFunctionBegin; 409 (*fact)->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqDense; 410 (*fact)->ops->solve = MatSolve_SeqDense; 411 (*fact)->ops->solvetranspose = MatSolveTranspose_SeqDense; 412 (*fact)->ops->solveadd = MatSolveAdd_SeqDense; 413 (*fact)->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense; 414 (*fact)->factor = MAT_FACTOR_CHOLESKY; 415 PetscFunctionReturn(0); 416 } 417 418 #undef __FUNCT__ 419 #define __FUNCT__ "MatLUFactorSymbolic_SeqDense" 420 PetscErrorCode MatLUFactorSymbolic_SeqDense(Mat A,IS row,IS col,MatFactorInfo *info,Mat *fact) 421 { 422 PetscFunctionBegin; 423 (*fact)->ops->lufactornumeric = MatLUFactorNumeric_SeqDense; 424 (*fact)->ops->solve = MatSolve_SeqDense; 425 (*fact)->ops->solvetranspose = MatSolveTranspose_SeqDense; 426 (*fact)->ops->solveadd = MatSolveAdd_SeqDense; 427 (*fact)->ops->solvetransposeadd = MatSolveTransposeAdd_SeqDense; 428 (*fact)->factor = MAT_FACTOR_LU; 429 PetscFunctionReturn(0); 430 } 431 432 #undef __FUNCT__ 433 #define __FUNCT__ "MatGetFactor_seqdense_petsc" 434 PetscErrorCode MatGetFactor_seqdense_petsc(Mat A,MatFactorType ftype,Mat *fact) 435 { 436 PetscErrorCode ierr; 437 438 PetscFunctionBegin; 439 ierr = MatCreate(((PetscObject)A)->comm,fact);CHKERRQ(ierr); 440 ierr = MatSetSizes(*fact,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 441 ierr = MatSetType(*fact,((PetscObject)A)->type_name);CHKERRQ(ierr); 442 if (ftype == MAT_FACTOR_LU){ 443 (*fact)->ops->lufactorsymbolic = MatLUFactorSymbolic_SeqDense; 444 } else { 445 (*fact)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqDense; 446 } 447 (*fact)->factor = ftype; 448 PetscFunctionReturn(0); 449 } 450 451 /* ------------------------------------------------------------------*/ 452 #undef __FUNCT__ 453 #define __FUNCT__ "MatRelax_SeqDense" 454 PetscErrorCode MatRelax_SeqDense(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal shift,PetscInt its,PetscInt lits,Vec xx) 455 { 456 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 457 PetscScalar *x,*b,*v = mat->v,zero = 0.0,xt; 458 PetscErrorCode ierr; 459 PetscInt m = A->rmap->n,i; 460 #if !defined(PETSC_USE_COMPLEX) 461 PetscBLASInt o = 1,bm = PetscBLASIntCast(m); 462 #endif 463 464 PetscFunctionBegin; 465 if (flag & SOR_ZERO_INITIAL_GUESS) { 466 /* this is a hack fix, should have another version without the second BLASdot */ 467 ierr = VecSet(xx,zero);CHKERRQ(ierr); 468 } 469 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 470 ierr = VecGetArray(bb,&b);CHKERRQ(ierr); 471 its = its*lits; 472 if (its <= 0) SETERRQ2(PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits); 473 while (its--) { 474 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 475 for (i=0; i<m; i++) { 476 #if defined(PETSC_USE_COMPLEX) 477 /* cannot use BLAS dot for complex because compiler/linker is 478 not happy about returning a double complex */ 479 PetscInt _i; 480 PetscScalar sum = b[i]; 481 for (_i=0; _i<m; _i++) { 482 sum -= PetscConj(v[i+_i*m])*x[_i]; 483 } 484 xt = sum; 485 #else 486 xt = b[i] - BLASdot_(&bm,v+i,&bm,x,&o); 487 #endif 488 x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift); 489 } 490 } 491 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 492 for (i=m-1; i>=0; i--) { 493 #if defined(PETSC_USE_COMPLEX) 494 /* cannot use BLAS dot for complex because compiler/linker is 495 not happy about returning a double complex */ 496 PetscInt _i; 497 PetscScalar sum = b[i]; 498 for (_i=0; _i<m; _i++) { 499 sum -= PetscConj(v[i+_i*m])*x[_i]; 500 } 501 xt = sum; 502 #else 503 xt = b[i] - BLASdot_(&bm,v+i,&bm,x,&o); 504 #endif 505 x[i] = (1. - omega)*x[i] + omega*(xt+v[i + i*m]*x[i])/(v[i + i*m]+shift); 506 } 507 } 508 } 509 ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr); 510 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 511 PetscFunctionReturn(0); 512 } 513 514 /* -----------------------------------------------------------------*/ 515 #undef __FUNCT__ 516 #define __FUNCT__ "MatMultTranspose_SeqDense" 517 PetscErrorCode MatMultTranspose_SeqDense(Mat A,Vec xx,Vec yy) 518 { 519 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 520 PetscScalar *v = mat->v,*x,*y; 521 PetscErrorCode ierr; 522 PetscBLASInt m, n,_One=1; 523 PetscScalar _DOne=1.0,_DZero=0.0; 524 525 PetscFunctionBegin; 526 m = PetscBLASIntCast(A->rmap->n); 527 n = PetscBLASIntCast(A->cmap->n); 528 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 529 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 530 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 531 BLASgemv_("T",&m,&n,&_DOne,v,&mat->lda,x,&_One,&_DZero,y,&_One); 532 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 533 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 534 ierr = PetscLogFlops(2*A->rmap->n*A->cmap->n - A->cmap->n);CHKERRQ(ierr); 535 PetscFunctionReturn(0); 536 } 537 538 #undef __FUNCT__ 539 #define __FUNCT__ "MatMult_SeqDense" 540 PetscErrorCode MatMult_SeqDense(Mat A,Vec xx,Vec yy) 541 { 542 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 543 PetscScalar *v = mat->v,*x,*y,_DOne=1.0,_DZero=0.0; 544 PetscErrorCode ierr; 545 PetscBLASInt m, n, _One=1; 546 547 PetscFunctionBegin; 548 m = PetscBLASIntCast(A->rmap->n); 549 n = PetscBLASIntCast(A->cmap->n); 550 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 551 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 552 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 553 BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DZero,y,&_One); 554 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 555 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 556 ierr = PetscLogFlops(2*A->rmap->n*A->cmap->n - A->rmap->n);CHKERRQ(ierr); 557 PetscFunctionReturn(0); 558 } 559 560 #undef __FUNCT__ 561 #define __FUNCT__ "MatMultAdd_SeqDense" 562 PetscErrorCode MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 563 { 564 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 565 PetscScalar *v = mat->v,*x,*y,_DOne=1.0; 566 PetscErrorCode ierr; 567 PetscBLASInt m, n, _One=1; 568 569 PetscFunctionBegin; 570 m = PetscBLASIntCast(A->rmap->n); 571 n = PetscBLASIntCast(A->cmap->n); 572 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 573 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 574 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 575 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 576 BLASgemv_("N",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One); 577 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 578 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 579 ierr = PetscLogFlops(2*A->rmap->n*A->cmap->n);CHKERRQ(ierr); 580 PetscFunctionReturn(0); 581 } 582 583 #undef __FUNCT__ 584 #define __FUNCT__ "MatMultTransposeAdd_SeqDense" 585 PetscErrorCode MatMultTransposeAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 586 { 587 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 588 PetscScalar *v = mat->v,*x,*y; 589 PetscErrorCode ierr; 590 PetscBLASInt m, n, _One=1; 591 PetscScalar _DOne=1.0; 592 593 PetscFunctionBegin; 594 m = PetscBLASIntCast(A->rmap->n); 595 n = PetscBLASIntCast(A->cmap->n); 596 if (!A->rmap->n || !A->cmap->n) PetscFunctionReturn(0); 597 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 598 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 599 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 600 BLASgemv_("T",&m,&n,&_DOne,v,&(mat->lda),x,&_One,&_DOne,y,&_One); 601 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 602 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 603 ierr = PetscLogFlops(2*A->rmap->n*A->cmap->n);CHKERRQ(ierr); 604 PetscFunctionReturn(0); 605 } 606 607 /* -----------------------------------------------------------------*/ 608 #undef __FUNCT__ 609 #define __FUNCT__ "MatGetRow_SeqDense" 610 PetscErrorCode MatGetRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals) 611 { 612 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 613 PetscScalar *v; 614 PetscErrorCode ierr; 615 PetscInt i; 616 617 PetscFunctionBegin; 618 *ncols = A->cmap->n; 619 if (cols) { 620 ierr = PetscMalloc((A->cmap->n+1)*sizeof(PetscInt),cols);CHKERRQ(ierr); 621 for (i=0; i<A->cmap->n; i++) (*cols)[i] = i; 622 } 623 if (vals) { 624 ierr = PetscMalloc((A->cmap->n+1)*sizeof(PetscScalar),vals);CHKERRQ(ierr); 625 v = mat->v + row; 626 for (i=0; i<A->cmap->n; i++) {(*vals)[i] = *v; v += mat->lda;} 627 } 628 PetscFunctionReturn(0); 629 } 630 631 #undef __FUNCT__ 632 #define __FUNCT__ "MatRestoreRow_SeqDense" 633 PetscErrorCode MatRestoreRow_SeqDense(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals) 634 { 635 PetscErrorCode ierr; 636 PetscFunctionBegin; 637 if (cols) {ierr = PetscFree(*cols);CHKERRQ(ierr);} 638 if (vals) {ierr = PetscFree(*vals);CHKERRQ(ierr); } 639 PetscFunctionReturn(0); 640 } 641 /* ----------------------------------------------------------------*/ 642 #undef __FUNCT__ 643 #define __FUNCT__ "MatSetValues_SeqDense" 644 PetscErrorCode MatSetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],const PetscScalar v[],InsertMode addv) 645 { 646 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 647 PetscInt i,j,idx=0; 648 649 PetscFunctionBegin; 650 if (!mat->roworiented) { 651 if (addv == INSERT_VALUES) { 652 for (j=0; j<n; j++) { 653 if (indexn[j] < 0) {idx += m; continue;} 654 #if defined(PETSC_USE_DEBUG) 655 if (indexn[j] >= A->cmap->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap->n-1); 656 #endif 657 for (i=0; i<m; i++) { 658 if (indexm[i] < 0) {idx++; continue;} 659 #if defined(PETSC_USE_DEBUG) 660 if (indexm[i] >= A->rmap->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap->n-1); 661 #endif 662 mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++]; 663 } 664 } 665 } else { 666 for (j=0; j<n; j++) { 667 if (indexn[j] < 0) {idx += m; continue;} 668 #if defined(PETSC_USE_DEBUG) 669 if (indexn[j] >= A->cmap->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap->n-1); 670 #endif 671 for (i=0; i<m; i++) { 672 if (indexm[i] < 0) {idx++; continue;} 673 #if defined(PETSC_USE_DEBUG) 674 if (indexm[i] >= A->rmap->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap->n-1); 675 #endif 676 mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++]; 677 } 678 } 679 } 680 } else { 681 if (addv == INSERT_VALUES) { 682 for (i=0; i<m; i++) { 683 if (indexm[i] < 0) { idx += n; continue;} 684 #if defined(PETSC_USE_DEBUG) 685 if (indexm[i] >= A->rmap->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap->n-1); 686 #endif 687 for (j=0; j<n; j++) { 688 if (indexn[j] < 0) { idx++; continue;} 689 #if defined(PETSC_USE_DEBUG) 690 if (indexn[j] >= A->cmap->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap->n-1); 691 #endif 692 mat->v[indexn[j]*mat->lda + indexm[i]] = v[idx++]; 693 } 694 } 695 } else { 696 for (i=0; i<m; i++) { 697 if (indexm[i] < 0) { idx += n; continue;} 698 #if defined(PETSC_USE_DEBUG) 699 if (indexm[i] >= A->rmap->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",indexm[i],A->rmap->n-1); 700 #endif 701 for (j=0; j<n; j++) { 702 if (indexn[j] < 0) { idx++; continue;} 703 #if defined(PETSC_USE_DEBUG) 704 if (indexn[j] >= A->cmap->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",indexn[j],A->cmap->n-1); 705 #endif 706 mat->v[indexn[j]*mat->lda + indexm[i]] += v[idx++]; 707 } 708 } 709 } 710 } 711 PetscFunctionReturn(0); 712 } 713 714 #undef __FUNCT__ 715 #define __FUNCT__ "MatGetValues_SeqDense" 716 PetscErrorCode MatGetValues_SeqDense(Mat A,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],PetscScalar v[]) 717 { 718 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 719 PetscInt i,j; 720 721 PetscFunctionBegin; 722 /* row-oriented output */ 723 for (i=0; i<m; i++) { 724 if (indexm[i] < 0) {v += n;continue;} 725 if (indexm[i] >= A->rmap->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row %D requested larger than number rows %D",indexm[i],A->rmap->n); 726 for (j=0; j<n; j++) { 727 if (indexn[j] < 0) {v++; continue;} 728 if (indexn[j] >= A->cmap->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column %D requested larger than number columns %D",indexn[j],A->cmap->n); 729 *v++ = mat->v[indexn[j]*mat->lda + indexm[i]]; 730 } 731 } 732 PetscFunctionReturn(0); 733 } 734 735 /* -----------------------------------------------------------------*/ 736 737 #include "petscsys.h" 738 739 #undef __FUNCT__ 740 #define __FUNCT__ "MatLoad_SeqDense" 741 PetscErrorCode MatLoad_SeqDense(PetscViewer viewer, const MatType type,Mat *A) 742 { 743 Mat_SeqDense *a; 744 Mat B; 745 PetscErrorCode ierr; 746 PetscInt *scols,i,j,nz,header[4]; 747 int fd; 748 PetscMPIInt size; 749 PetscInt *rowlengths = 0,M,N,*cols; 750 PetscScalar *vals,*svals,*v,*w; 751 MPI_Comm comm = ((PetscObject)viewer)->comm; 752 753 PetscFunctionBegin; 754 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 755 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"view must have one processor"); 756 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 757 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 758 if (header[0] != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Not matrix object"); 759 M = header[1]; N = header[2]; nz = header[3]; 760 761 if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense */ 762 ierr = MatCreate(comm,A);CHKERRQ(ierr); 763 ierr = MatSetSizes(*A,M,N,M,N);CHKERRQ(ierr); 764 ierr = MatSetType(*A,type);CHKERRQ(ierr); 765 ierr = MatSeqDenseSetPreallocation(*A,PETSC_NULL);CHKERRQ(ierr); 766 B = *A; 767 a = (Mat_SeqDense*)B->data; 768 v = a->v; 769 /* Allocate some temp space to read in the values and then flip them 770 from row major to column major */ 771 ierr = PetscMalloc((M*N > 0 ? M*N : 1)*sizeof(PetscScalar),&w);CHKERRQ(ierr); 772 /* read in nonzero values */ 773 ierr = PetscBinaryRead(fd,w,M*N,PETSC_SCALAR);CHKERRQ(ierr); 774 /* now flip the values and store them in the matrix*/ 775 for (j=0; j<N; j++) { 776 for (i=0; i<M; i++) { 777 *v++ =w[i*N+j]; 778 } 779 } 780 ierr = PetscFree(w);CHKERRQ(ierr); 781 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 782 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 783 } else { 784 /* read row lengths */ 785 ierr = PetscMalloc((M+1)*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 786 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 787 788 /* create our matrix */ 789 ierr = MatCreate(comm,A);CHKERRQ(ierr); 790 ierr = MatSetSizes(*A,M,N,M,N);CHKERRQ(ierr); 791 ierr = MatSetType(*A,type);CHKERRQ(ierr); 792 ierr = MatSeqDenseSetPreallocation(*A,PETSC_NULL);CHKERRQ(ierr); 793 B = *A; 794 a = (Mat_SeqDense*)B->data; 795 v = a->v; 796 797 /* read column indices and nonzeros */ 798 ierr = PetscMalloc((nz+1)*sizeof(PetscInt),&scols);CHKERRQ(ierr); 799 cols = scols; 800 ierr = PetscBinaryRead(fd,cols,nz,PETSC_INT);CHKERRQ(ierr); 801 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&svals);CHKERRQ(ierr); 802 vals = svals; 803 ierr = PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);CHKERRQ(ierr); 804 805 /* insert into matrix */ 806 for (i=0; i<M; i++) { 807 for (j=0; j<rowlengths[i]; j++) v[i+M*scols[j]] = svals[j]; 808 svals += rowlengths[i]; scols += rowlengths[i]; 809 } 810 ierr = PetscFree(vals);CHKERRQ(ierr); 811 ierr = PetscFree(cols);CHKERRQ(ierr); 812 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 813 814 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 815 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 816 } 817 PetscFunctionReturn(0); 818 } 819 820 #include "petscsys.h" 821 822 #undef __FUNCT__ 823 #define __FUNCT__ "MatView_SeqDense_ASCII" 824 static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer) 825 { 826 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 827 PetscErrorCode ierr; 828 PetscInt i,j; 829 const char *name; 830 PetscScalar *v; 831 PetscViewerFormat format; 832 #if defined(PETSC_USE_COMPLEX) 833 PetscTruth allreal = PETSC_TRUE; 834 #endif 835 836 PetscFunctionBegin; 837 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 838 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 839 if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 840 PetscFunctionReturn(0); /* do nothing for now */ 841 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 842 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 843 for (i=0; i<A->rmap->n; i++) { 844 v = a->v + i; 845 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 846 for (j=0; j<A->cmap->n; j++) { 847 #if defined(PETSC_USE_COMPLEX) 848 if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) { 849 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i) ",j,PetscRealPart(*v),PetscImaginaryPart(*v));CHKERRQ(ierr); 850 } else if (PetscRealPart(*v)) { 851 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",j,PetscRealPart(*v));CHKERRQ(ierr); 852 } 853 #else 854 if (*v) { 855 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",j,*v);CHKERRQ(ierr); 856 } 857 #endif 858 v += a->lda; 859 } 860 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 861 } 862 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 863 } else { 864 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 865 #if defined(PETSC_USE_COMPLEX) 866 /* determine if matrix has all real values */ 867 v = a->v; 868 for (i=0; i<A->rmap->n*A->cmap->n; i++) { 869 if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break ;} 870 } 871 #endif 872 if (format == PETSC_VIEWER_ASCII_MATLAB) { 873 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 874 ierr = PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",A->rmap->n,A->cmap->n);CHKERRQ(ierr); 875 ierr = PetscViewerASCIIPrintf(viewer,"%s = zeros(%D,%D);\n",name,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 876 ierr = PetscViewerASCIIPrintf(viewer,"%s = [\n",name);CHKERRQ(ierr); 877 } 878 879 for (i=0; i<A->rmap->n; i++) { 880 v = a->v + i; 881 for (j=0; j<A->cmap->n; j++) { 882 #if defined(PETSC_USE_COMPLEX) 883 if (allreal) { 884 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e ",PetscRealPart(*v));CHKERRQ(ierr); 885 } else { 886 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16e i ",PetscRealPart(*v),PetscImaginaryPart(*v));CHKERRQ(ierr); 887 } 888 #else 889 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e ",*v);CHKERRQ(ierr); 890 #endif 891 v += a->lda; 892 } 893 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 894 } 895 if (format == PETSC_VIEWER_ASCII_MATLAB) { 896 ierr = PetscViewerASCIIPrintf(viewer,"];\n");CHKERRQ(ierr); 897 } 898 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 899 } 900 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 901 PetscFunctionReturn(0); 902 } 903 904 #undef __FUNCT__ 905 #define __FUNCT__ "MatView_SeqDense_Binary" 906 static PetscErrorCode MatView_SeqDense_Binary(Mat A,PetscViewer viewer) 907 { 908 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 909 PetscErrorCode ierr; 910 int fd; 911 PetscInt ict,j,n = A->cmap->n,m = A->rmap->n,i,*col_lens,nz = m*n; 912 PetscScalar *v,*anonz,*vals; 913 PetscViewerFormat format; 914 915 PetscFunctionBegin; 916 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 917 918 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 919 if (format == PETSC_VIEWER_BINARY_NATIVE) { 920 /* store the matrix as a dense matrix */ 921 ierr = PetscMalloc(4*sizeof(PetscInt),&col_lens);CHKERRQ(ierr); 922 col_lens[0] = MAT_FILE_COOKIE; 923 col_lens[1] = m; 924 col_lens[2] = n; 925 col_lens[3] = MATRIX_BINARY_FORMAT_DENSE; 926 ierr = PetscBinaryWrite(fd,col_lens,4,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 927 ierr = PetscFree(col_lens);CHKERRQ(ierr); 928 929 /* write out matrix, by rows */ 930 ierr = PetscMalloc((m*n+1)*sizeof(PetscScalar),&vals);CHKERRQ(ierr); 931 v = a->v; 932 for (j=0; j<n; j++) { 933 for (i=0; i<m; i++) { 934 vals[j + i*n] = *v++; 935 } 936 } 937 ierr = PetscBinaryWrite(fd,vals,n*m,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 938 ierr = PetscFree(vals);CHKERRQ(ierr); 939 } else { 940 ierr = PetscMalloc((4+nz)*sizeof(PetscInt),&col_lens);CHKERRQ(ierr); 941 col_lens[0] = MAT_FILE_COOKIE; 942 col_lens[1] = m; 943 col_lens[2] = n; 944 col_lens[3] = nz; 945 946 /* store lengths of each row and write (including header) to file */ 947 for (i=0; i<m; i++) col_lens[4+i] = n; 948 ierr = PetscBinaryWrite(fd,col_lens,4+m,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 949 950 /* Possibly should write in smaller increments, not whole matrix at once? */ 951 /* store column indices (zero start index) */ 952 ict = 0; 953 for (i=0; i<m; i++) { 954 for (j=0; j<n; j++) col_lens[ict++] = j; 955 } 956 ierr = PetscBinaryWrite(fd,col_lens,nz,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 957 ierr = PetscFree(col_lens);CHKERRQ(ierr); 958 959 /* store nonzero values */ 960 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&anonz);CHKERRQ(ierr); 961 ict = 0; 962 for (i=0; i<m; i++) { 963 v = a->v + i; 964 for (j=0; j<n; j++) { 965 anonz[ict++] = *v; v += a->lda; 966 } 967 } 968 ierr = PetscBinaryWrite(fd,anonz,nz,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 969 ierr = PetscFree(anonz);CHKERRQ(ierr); 970 } 971 PetscFunctionReturn(0); 972 } 973 974 #undef __FUNCT__ 975 #define __FUNCT__ "MatView_SeqDense_Draw_Zoom" 976 PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa) 977 { 978 Mat A = (Mat) Aa; 979 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 980 PetscErrorCode ierr; 981 PetscInt m = A->rmap->n,n = A->cmap->n,color,i,j; 982 PetscScalar *v = a->v; 983 PetscViewer viewer; 984 PetscDraw popup; 985 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r,scale,maxv = 0.0; 986 PetscViewerFormat format; 987 988 PetscFunctionBegin; 989 990 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 991 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 992 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 993 994 /* Loop over matrix elements drawing boxes */ 995 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 996 /* Blue for negative and Red for positive */ 997 color = PETSC_DRAW_BLUE; 998 for(j = 0; j < n; j++) { 999 x_l = j; 1000 x_r = x_l + 1.0; 1001 for(i = 0; i < m; i++) { 1002 y_l = m - i - 1.0; 1003 y_r = y_l + 1.0; 1004 #if defined(PETSC_USE_COMPLEX) 1005 if (PetscRealPart(v[j*m+i]) > 0.) { 1006 color = PETSC_DRAW_RED; 1007 } else if (PetscRealPart(v[j*m+i]) < 0.) { 1008 color = PETSC_DRAW_BLUE; 1009 } else { 1010 continue; 1011 } 1012 #else 1013 if (v[j*m+i] > 0.) { 1014 color = PETSC_DRAW_RED; 1015 } else if (v[j*m+i] < 0.) { 1016 color = PETSC_DRAW_BLUE; 1017 } else { 1018 continue; 1019 } 1020 #endif 1021 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 1022 } 1023 } 1024 } else { 1025 /* use contour shading to indicate magnitude of values */ 1026 /* first determine max of all nonzero values */ 1027 for(i = 0; i < m*n; i++) { 1028 if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]); 1029 } 1030 scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/maxv; 1031 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 1032 if (popup) {ierr = PetscDrawScalePopup(popup,0.0,maxv);CHKERRQ(ierr);} 1033 for(j = 0; j < n; j++) { 1034 x_l = j; 1035 x_r = x_l + 1.0; 1036 for(i = 0; i < m; i++) { 1037 y_l = m - i - 1.0; 1038 y_r = y_l + 1.0; 1039 color = PETSC_DRAW_BASIC_COLORS + (int)(scale*PetscAbsScalar(v[j*m+i])); 1040 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 1041 } 1042 } 1043 } 1044 PetscFunctionReturn(0); 1045 } 1046 1047 #undef __FUNCT__ 1048 #define __FUNCT__ "MatView_SeqDense_Draw" 1049 PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer) 1050 { 1051 PetscDraw draw; 1052 PetscTruth isnull; 1053 PetscReal xr,yr,xl,yl,h,w; 1054 PetscErrorCode ierr; 1055 1056 PetscFunctionBegin; 1057 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 1058 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 1059 if (isnull) PetscFunctionReturn(0); 1060 1061 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 1062 xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0; 1063 xr += w; yr += h; xl = -w; yl = -h; 1064 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 1065 ierr = PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);CHKERRQ(ierr); 1066 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);CHKERRQ(ierr); 1067 PetscFunctionReturn(0); 1068 } 1069 1070 #undef __FUNCT__ 1071 #define __FUNCT__ "MatView_SeqDense" 1072 PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer) 1073 { 1074 PetscErrorCode ierr; 1075 PetscTruth iascii,isbinary,isdraw; 1076 1077 PetscFunctionBegin; 1078 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 1079 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_BINARY,&isbinary);CHKERRQ(ierr); 1080 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_DRAW,&isdraw);CHKERRQ(ierr); 1081 1082 if (iascii) { 1083 ierr = MatView_SeqDense_ASCII(A,viewer);CHKERRQ(ierr); 1084 } else if (isbinary) { 1085 ierr = MatView_SeqDense_Binary(A,viewer);CHKERRQ(ierr); 1086 } else if (isdraw) { 1087 ierr = MatView_SeqDense_Draw(A,viewer);CHKERRQ(ierr); 1088 } else { 1089 SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported by dense matrix",((PetscObject)viewer)->type_name); 1090 } 1091 PetscFunctionReturn(0); 1092 } 1093 1094 #undef __FUNCT__ 1095 #define __FUNCT__ "MatDestroy_SeqDense" 1096 PetscErrorCode MatDestroy_SeqDense(Mat mat) 1097 { 1098 Mat_SeqDense *l = (Mat_SeqDense*)mat->data; 1099 PetscErrorCode ierr; 1100 1101 PetscFunctionBegin; 1102 #if defined(PETSC_USE_LOG) 1103 PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->rmap->n,mat->cmap->n); 1104 #endif 1105 ierr = PetscFree(l->pivots);CHKERRQ(ierr); 1106 if (!l->user_alloc) {ierr = PetscFree(l->v);CHKERRQ(ierr);} 1107 ierr = PetscFree(l);CHKERRQ(ierr); 1108 1109 ierr = PetscObjectChangeTypeName((PetscObject)mat,0);CHKERRQ(ierr); 1110 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatSeqDenseSetPreallocation_C","",PETSC_NULL);CHKERRQ(ierr); 1111 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatMatMult_seqaij_seqdense_C","",PETSC_NULL);CHKERRQ(ierr); 1112 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatMatMultSymbolic_seqaij_seqdense_C","",PETSC_NULL);CHKERRQ(ierr); 1113 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatMatMultNumeric_seqaij_seqdense_C","",PETSC_NULL);CHKERRQ(ierr); 1114 PetscFunctionReturn(0); 1115 } 1116 1117 #undef __FUNCT__ 1118 #define __FUNCT__ "MatTranspose_SeqDense" 1119 PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout) 1120 { 1121 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1122 PetscErrorCode ierr; 1123 PetscInt k,j,m,n,M; 1124 PetscScalar *v,tmp; 1125 1126 PetscFunctionBegin; 1127 v = mat->v; m = A->rmap->n; M = mat->lda; n = A->cmap->n; 1128 if (reuse == MAT_REUSE_MATRIX && *matout == A) { /* in place transpose */ 1129 if (m != n) { 1130 SETERRQ(PETSC_ERR_SUP,"Can not transpose non-square matrix in place"); 1131 } else { 1132 for (j=0; j<m; j++) { 1133 for (k=0; k<j; k++) { 1134 tmp = v[j + k*M]; 1135 v[j + k*M] = v[k + j*M]; 1136 v[k + j*M] = tmp; 1137 } 1138 } 1139 } 1140 } else { /* out-of-place transpose */ 1141 Mat tmat; 1142 Mat_SeqDense *tmatd; 1143 PetscScalar *v2; 1144 1145 if (reuse == MAT_INITIAL_MATRIX) { 1146 ierr = MatCreate(((PetscObject)A)->comm,&tmat);CHKERRQ(ierr); 1147 ierr = MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);CHKERRQ(ierr); 1148 ierr = MatSetType(tmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1149 ierr = MatSeqDenseSetPreallocation(tmat,PETSC_NULL);CHKERRQ(ierr); 1150 } else { 1151 tmat = *matout; 1152 } 1153 tmatd = (Mat_SeqDense*)tmat->data; 1154 v = mat->v; v2 = tmatd->v; 1155 for (j=0; j<n; j++) { 1156 for (k=0; k<m; k++) v2[j + k*n] = v[k + j*M]; 1157 } 1158 ierr = MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1159 ierr = MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1160 *matout = tmat; 1161 } 1162 PetscFunctionReturn(0); 1163 } 1164 1165 #undef __FUNCT__ 1166 #define __FUNCT__ "MatEqual_SeqDense" 1167 PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscTruth *flg) 1168 { 1169 Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data; 1170 Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data; 1171 PetscInt i,j; 1172 PetscScalar *v1 = mat1->v,*v2 = mat2->v; 1173 1174 PetscFunctionBegin; 1175 if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1176 if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1177 for (i=0; i<A1->rmap->n; i++) { 1178 v1 = mat1->v+i; v2 = mat2->v+i; 1179 for (j=0; j<A1->cmap->n; j++) { 1180 if (*v1 != *v2) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1181 v1 += mat1->lda; v2 += mat2->lda; 1182 } 1183 } 1184 *flg = PETSC_TRUE; 1185 PetscFunctionReturn(0); 1186 } 1187 1188 #undef __FUNCT__ 1189 #define __FUNCT__ "MatGetDiagonal_SeqDense" 1190 PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v) 1191 { 1192 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1193 PetscErrorCode ierr; 1194 PetscInt i,n,len; 1195 PetscScalar *x,zero = 0.0; 1196 1197 PetscFunctionBegin; 1198 ierr = VecSet(v,zero);CHKERRQ(ierr); 1199 ierr = VecGetSize(v,&n);CHKERRQ(ierr); 1200 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1201 len = PetscMin(A->rmap->n,A->cmap->n); 1202 if (n != A->rmap->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec"); 1203 for (i=0; i<len; i++) { 1204 x[i] = mat->v[i*mat->lda + i]; 1205 } 1206 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1207 PetscFunctionReturn(0); 1208 } 1209 1210 #undef __FUNCT__ 1211 #define __FUNCT__ "MatDiagonalScale_SeqDense" 1212 PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr) 1213 { 1214 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1215 PetscScalar *l,*r,x,*v; 1216 PetscErrorCode ierr; 1217 PetscInt i,j,m = A->rmap->n,n = A->cmap->n; 1218 1219 PetscFunctionBegin; 1220 if (ll) { 1221 ierr = VecGetSize(ll,&m);CHKERRQ(ierr); 1222 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1223 if (m != A->rmap->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size"); 1224 for (i=0; i<m; i++) { 1225 x = l[i]; 1226 v = mat->v + i; 1227 for (j=0; j<n; j++) { (*v) *= x; v+= m;} 1228 } 1229 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1230 ierr = PetscLogFlops(n*m);CHKERRQ(ierr); 1231 } 1232 if (rr) { 1233 ierr = VecGetSize(rr,&n);CHKERRQ(ierr); 1234 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1235 if (n != A->cmap->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size"); 1236 for (i=0; i<n; i++) { 1237 x = r[i]; 1238 v = mat->v + i*m; 1239 for (j=0; j<m; j++) { (*v++) *= x;} 1240 } 1241 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1242 ierr = PetscLogFlops(n*m);CHKERRQ(ierr); 1243 } 1244 PetscFunctionReturn(0); 1245 } 1246 1247 #undef __FUNCT__ 1248 #define __FUNCT__ "MatNorm_SeqDense" 1249 PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm) 1250 { 1251 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1252 PetscScalar *v = mat->v; 1253 PetscReal sum = 0.0; 1254 PetscInt lda=mat->lda,m=A->rmap->n,i,j; 1255 PetscErrorCode ierr; 1256 1257 PetscFunctionBegin; 1258 if (type == NORM_FROBENIUS) { 1259 if (lda>m) { 1260 for (j=0; j<A->cmap->n; j++) { 1261 v = mat->v+j*lda; 1262 for (i=0; i<m; i++) { 1263 #if defined(PETSC_USE_COMPLEX) 1264 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1265 #else 1266 sum += (*v)*(*v); v++; 1267 #endif 1268 } 1269 } 1270 } else { 1271 for (i=0; i<A->cmap->n*A->rmap->n; i++) { 1272 #if defined(PETSC_USE_COMPLEX) 1273 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1274 #else 1275 sum += (*v)*(*v); v++; 1276 #endif 1277 } 1278 } 1279 *nrm = sqrt(sum); 1280 ierr = PetscLogFlops(2*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1281 } else if (type == NORM_1) { 1282 *nrm = 0.0; 1283 for (j=0; j<A->cmap->n; j++) { 1284 v = mat->v + j*mat->lda; 1285 sum = 0.0; 1286 for (i=0; i<A->rmap->n; i++) { 1287 sum += PetscAbsScalar(*v); v++; 1288 } 1289 if (sum > *nrm) *nrm = sum; 1290 } 1291 ierr = PetscLogFlops(A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1292 } else if (type == NORM_INFINITY) { 1293 *nrm = 0.0; 1294 for (j=0; j<A->rmap->n; j++) { 1295 v = mat->v + j; 1296 sum = 0.0; 1297 for (i=0; i<A->cmap->n; i++) { 1298 sum += PetscAbsScalar(*v); v += mat->lda; 1299 } 1300 if (sum > *nrm) *nrm = sum; 1301 } 1302 ierr = PetscLogFlops(A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1303 } else { 1304 SETERRQ(PETSC_ERR_SUP,"No two norm"); 1305 } 1306 PetscFunctionReturn(0); 1307 } 1308 1309 #undef __FUNCT__ 1310 #define __FUNCT__ "MatSetOption_SeqDense" 1311 PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscTruth flg) 1312 { 1313 Mat_SeqDense *aij = (Mat_SeqDense*)A->data; 1314 PetscErrorCode ierr; 1315 1316 PetscFunctionBegin; 1317 switch (op) { 1318 case MAT_ROW_ORIENTED: 1319 aij->roworiented = flg; 1320 break; 1321 case MAT_NEW_NONZERO_LOCATIONS: 1322 case MAT_NEW_NONZERO_LOCATION_ERR: 1323 case MAT_NEW_NONZERO_ALLOCATION_ERR: 1324 case MAT_NEW_DIAGONALS: 1325 case MAT_IGNORE_OFF_PROC_ENTRIES: 1326 case MAT_USE_HASH_TABLE: 1327 case MAT_SYMMETRIC: 1328 case MAT_STRUCTURALLY_SYMMETRIC: 1329 case MAT_HERMITIAN: 1330 case MAT_SYMMETRY_ETERNAL: 1331 case MAT_IGNORE_LOWER_TRIANGULAR: 1332 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 1333 break; 1334 default: 1335 SETERRQ1(PETSC_ERR_SUP,"unknown option %s",MatOptions[op]); 1336 } 1337 PetscFunctionReturn(0); 1338 } 1339 1340 #undef __FUNCT__ 1341 #define __FUNCT__ "MatZeroEntries_SeqDense" 1342 PetscErrorCode MatZeroEntries_SeqDense(Mat A) 1343 { 1344 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1345 PetscErrorCode ierr; 1346 PetscInt lda=l->lda,m=A->rmap->n,j; 1347 1348 PetscFunctionBegin; 1349 if (lda>m) { 1350 for (j=0; j<A->cmap->n; j++) { 1351 ierr = PetscMemzero(l->v+j*lda,m*sizeof(PetscScalar));CHKERRQ(ierr); 1352 } 1353 } else { 1354 ierr = PetscMemzero(l->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1355 } 1356 PetscFunctionReturn(0); 1357 } 1358 1359 #undef __FUNCT__ 1360 #define __FUNCT__ "MatZeroRows_SeqDense" 1361 PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag) 1362 { 1363 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1364 PetscInt n = A->cmap->n,i,j; 1365 PetscScalar *slot; 1366 1367 PetscFunctionBegin; 1368 for (i=0; i<N; i++) { 1369 slot = l->v + rows[i]; 1370 for (j=0; j<n; j++) { *slot = 0.0; slot += n;} 1371 } 1372 if (diag != 0.0) { 1373 for (i=0; i<N; i++) { 1374 slot = l->v + (n+1)*rows[i]; 1375 *slot = diag; 1376 } 1377 } 1378 PetscFunctionReturn(0); 1379 } 1380 1381 #undef __FUNCT__ 1382 #define __FUNCT__ "MatGetArray_SeqDense" 1383 PetscErrorCode MatGetArray_SeqDense(Mat A,PetscScalar *array[]) 1384 { 1385 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1386 1387 PetscFunctionBegin; 1388 if (mat->lda != A->rmap->n) SETERRQ(PETSC_ERR_SUP,"Cannot get array for Dense matrices with LDA different from number of rows"); 1389 *array = mat->v; 1390 PetscFunctionReturn(0); 1391 } 1392 1393 #undef __FUNCT__ 1394 #define __FUNCT__ "MatRestoreArray_SeqDense" 1395 PetscErrorCode MatRestoreArray_SeqDense(Mat A,PetscScalar *array[]) 1396 { 1397 PetscFunctionBegin; 1398 *array = 0; /* user cannot accidently use the array later */ 1399 PetscFunctionReturn(0); 1400 } 1401 1402 #undef __FUNCT__ 1403 #define __FUNCT__ "MatGetSubMatrix_SeqDense" 1404 static PetscErrorCode MatGetSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B) 1405 { 1406 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1407 PetscErrorCode ierr; 1408 PetscInt i,j,*irow,*icol,nrows,ncols; 1409 PetscScalar *av,*bv,*v = mat->v; 1410 Mat newmat; 1411 1412 PetscFunctionBegin; 1413 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1414 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1415 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1416 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1417 1418 /* Check submatrixcall */ 1419 if (scall == MAT_REUSE_MATRIX) { 1420 PetscInt n_cols,n_rows; 1421 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1422 if (n_rows != nrows || n_cols != ncols) { 1423 /* resize the result result matrix to match number of requested rows/columns */ 1424 ierr = MatSetSizes(*B,nrows,nrows,nrows,nrows);CHKERRQ(ierr); 1425 } 1426 newmat = *B; 1427 } else { 1428 /* Create and fill new matrix */ 1429 ierr = MatCreate(((PetscObject)A)->comm,&newmat);CHKERRQ(ierr); 1430 ierr = MatSetSizes(newmat,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1431 ierr = MatSetType(newmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1432 ierr = MatSeqDenseSetPreallocation(newmat,PETSC_NULL);CHKERRQ(ierr); 1433 } 1434 1435 /* Now extract the data pointers and do the copy,column at a time */ 1436 bv = ((Mat_SeqDense*)newmat->data)->v; 1437 1438 for (i=0; i<ncols; i++) { 1439 av = v + mat->lda*icol[i]; 1440 for (j=0; j<nrows; j++) { 1441 *bv++ = av[irow[j]]; 1442 } 1443 } 1444 1445 /* Assemble the matrices so that the correct flags are set */ 1446 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1447 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1448 1449 /* Free work space */ 1450 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1451 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1452 *B = newmat; 1453 PetscFunctionReturn(0); 1454 } 1455 1456 #undef __FUNCT__ 1457 #define __FUNCT__ "MatGetSubMatrices_SeqDense" 1458 PetscErrorCode MatGetSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1459 { 1460 PetscErrorCode ierr; 1461 PetscInt i; 1462 1463 PetscFunctionBegin; 1464 if (scall == MAT_INITIAL_MATRIX) { 1465 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1466 } 1467 1468 for (i=0; i<n; i++) { 1469 ierr = MatGetSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1470 } 1471 PetscFunctionReturn(0); 1472 } 1473 1474 #undef __FUNCT__ 1475 #define __FUNCT__ "MatAssemblyBegin_SeqDense" 1476 PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode) 1477 { 1478 PetscFunctionBegin; 1479 PetscFunctionReturn(0); 1480 } 1481 1482 #undef __FUNCT__ 1483 #define __FUNCT__ "MatAssemblyEnd_SeqDense" 1484 PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode) 1485 { 1486 PetscFunctionBegin; 1487 PetscFunctionReturn(0); 1488 } 1489 1490 #undef __FUNCT__ 1491 #define __FUNCT__ "MatCopy_SeqDense" 1492 PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str) 1493 { 1494 Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense *)B->data; 1495 PetscErrorCode ierr; 1496 PetscInt lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j; 1497 1498 PetscFunctionBegin; 1499 /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */ 1500 if (A->ops->copy != B->ops->copy) { 1501 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1502 PetscFunctionReturn(0); 1503 } 1504 if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_ERR_ARG_SIZ,"size(B) != size(A)"); 1505 if (lda1>m || lda2>m) { 1506 for (j=0; j<n; j++) { 1507 ierr = PetscMemcpy(b->v+j*lda2,a->v+j*lda1,m*sizeof(PetscScalar));CHKERRQ(ierr); 1508 } 1509 } else { 1510 ierr = PetscMemcpy(b->v,a->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1511 } 1512 PetscFunctionReturn(0); 1513 } 1514 1515 #undef __FUNCT__ 1516 #define __FUNCT__ "MatSetUpPreallocation_SeqDense" 1517 PetscErrorCode MatSetUpPreallocation_SeqDense(Mat A) 1518 { 1519 PetscErrorCode ierr; 1520 1521 PetscFunctionBegin; 1522 ierr = MatSeqDenseSetPreallocation(A,0);CHKERRQ(ierr); 1523 PetscFunctionReturn(0); 1524 } 1525 1526 #undef __FUNCT__ 1527 #define __FUNCT__ "MatSetSizes_SeqDense" 1528 PetscErrorCode MatSetSizes_SeqDense(Mat A,PetscInt m,PetscInt n,PetscInt M,PetscInt N) 1529 { 1530 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1531 PetscErrorCode ierr; 1532 PetscFunctionBegin; 1533 /* this will not be called before lda, Mmax, and Nmax have been set */ 1534 m = PetscMax(m,M); 1535 n = PetscMax(n,N); 1536 if (m > a->Mmax) SETERRQ2(PETSC_ERR_SUP,"Cannot yet resize number rows of dense matrix larger then its initial size %d, requested %d",a->lda,(int)m); 1537 if (n > a->Nmax) SETERRQ2(PETSC_ERR_SUP,"Cannot yet resize number columns of dense matrix larger then its initial size %d, requested %d",a->Nmax,(int)n); 1538 A->rmap->n = A->rmap->n = m; 1539 A->cmap->n = A->cmap->N = n; 1540 if (a->changelda) a->lda = m; 1541 ierr = PetscMemzero(a->v,m*n*sizeof(PetscScalar));CHKERRQ(ierr); 1542 PetscFunctionReturn(0); 1543 } 1544 1545 1546 /* ----------------------------------------------------------------*/ 1547 #undef __FUNCT__ 1548 #define __FUNCT__ "MatMatMult_SeqDense_SeqDense" 1549 PetscErrorCode MatMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 1550 { 1551 PetscErrorCode ierr; 1552 1553 PetscFunctionBegin; 1554 if (scall == MAT_INITIAL_MATRIX){ 1555 ierr = MatMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 1556 } 1557 ierr = MatMatMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 1558 PetscFunctionReturn(0); 1559 } 1560 1561 #undef __FUNCT__ 1562 #define __FUNCT__ "MatMatMultSymbolic_SeqDense_SeqDense" 1563 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 1564 { 1565 PetscErrorCode ierr; 1566 PetscInt m=A->rmap->n,n=B->cmap->n; 1567 Mat Cmat; 1568 1569 PetscFunctionBegin; 1570 if (A->cmap->n != B->rmap->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"A->cmap->n %d != B->rmap->n %d\n",A->cmap->n,B->rmap->n); 1571 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 1572 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 1573 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 1574 ierr = MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);CHKERRQ(ierr); 1575 Cmat->assembled = PETSC_TRUE; 1576 *C = Cmat; 1577 PetscFunctionReturn(0); 1578 } 1579 1580 #undef __FUNCT__ 1581 #define __FUNCT__ "MatMatMultNumeric_SeqDense_SeqDense" 1582 PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 1583 { 1584 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1585 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1586 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 1587 PetscBLASInt m,n,k; 1588 PetscScalar _DOne=1.0,_DZero=0.0; 1589 1590 PetscFunctionBegin; 1591 m = PetscBLASIntCast(A->rmap->n); 1592 n = PetscBLASIntCast(B->cmap->n); 1593 k = PetscBLASIntCast(A->cmap->n); 1594 BLASgemm_("N","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda); 1595 PetscFunctionReturn(0); 1596 } 1597 1598 #undef __FUNCT__ 1599 #define __FUNCT__ "MatMatMultTranspose_SeqDense_SeqDense" 1600 PetscErrorCode MatMatMultTranspose_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 1601 { 1602 PetscErrorCode ierr; 1603 1604 PetscFunctionBegin; 1605 if (scall == MAT_INITIAL_MATRIX){ 1606 ierr = MatMatMultTransposeSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 1607 } 1608 ierr = MatMatMultTransposeNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 1609 PetscFunctionReturn(0); 1610 } 1611 1612 #undef __FUNCT__ 1613 #define __FUNCT__ "MatMatMultTransposeSymbolic_SeqDense_SeqDense" 1614 PetscErrorCode MatMatMultTransposeSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 1615 { 1616 PetscErrorCode ierr; 1617 PetscInt m=A->cmap->n,n=B->cmap->n; 1618 Mat Cmat; 1619 1620 PetscFunctionBegin; 1621 if (A->rmap->n != B->rmap->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"A->rmap->n %d != B->rmap->n %d\n",A->rmap->n,B->rmap->n); 1622 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 1623 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 1624 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 1625 ierr = MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);CHKERRQ(ierr); 1626 Cmat->assembled = PETSC_TRUE; 1627 *C = Cmat; 1628 PetscFunctionReturn(0); 1629 } 1630 1631 #undef __FUNCT__ 1632 #define __FUNCT__ "MatMatMultTransposeNumeric_SeqDense_SeqDense" 1633 PetscErrorCode MatMatMultTransposeNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 1634 { 1635 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1636 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1637 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 1638 PetscBLASInt m,n,k; 1639 PetscScalar _DOne=1.0,_DZero=0.0; 1640 1641 PetscFunctionBegin; 1642 m = PetscBLASIntCast(A->cmap->n); 1643 n = PetscBLASIntCast(B->cmap->n); 1644 k = PetscBLASIntCast(A->rmap->n); 1645 /* 1646 Note the m and n arguments below are the number rows and columns of A', not A! 1647 */ 1648 BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda); 1649 PetscFunctionReturn(0); 1650 } 1651 1652 #undef __FUNCT__ 1653 #define __FUNCT__ "MatGetRowMax_SeqDense" 1654 PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[]) 1655 { 1656 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1657 PetscErrorCode ierr; 1658 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1659 PetscScalar *x; 1660 MatScalar *aa = a->v; 1661 1662 PetscFunctionBegin; 1663 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1664 1665 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1666 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1667 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1668 if (p != A->rmap->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1669 for (i=0; i<m; i++) { 1670 x[i] = aa[i]; if (idx) idx[i] = 0; 1671 for (j=1; j<n; j++){ 1672 if (PetscRealPart(x[i]) < PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;} 1673 } 1674 } 1675 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1676 PetscFunctionReturn(0); 1677 } 1678 1679 #undef __FUNCT__ 1680 #define __FUNCT__ "MatGetRowMaxAbs_SeqDense" 1681 PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[]) 1682 { 1683 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1684 PetscErrorCode ierr; 1685 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1686 PetscScalar *x; 1687 PetscReal atmp; 1688 MatScalar *aa = a->v; 1689 1690 PetscFunctionBegin; 1691 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1692 1693 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1694 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1695 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1696 if (p != A->rmap->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1697 for (i=0; i<m; i++) { 1698 x[i] = PetscAbsScalar(aa[i]); 1699 for (j=1; j<n; j++){ 1700 atmp = PetscAbsScalar(aa[i+m*j]); 1701 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;} 1702 } 1703 } 1704 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1705 PetscFunctionReturn(0); 1706 } 1707 1708 #undef __FUNCT__ 1709 #define __FUNCT__ "MatGetRowMin_SeqDense" 1710 PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[]) 1711 { 1712 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1713 PetscErrorCode ierr; 1714 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1715 PetscScalar *x; 1716 MatScalar *aa = a->v; 1717 1718 PetscFunctionBegin; 1719 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1720 1721 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1722 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1723 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1724 if (p != A->rmap->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1725 for (i=0; i<m; i++) { 1726 x[i] = aa[i]; if (idx) idx[i] = 0; 1727 for (j=1; j<n; j++){ 1728 if (PetscRealPart(x[i]) > PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;} 1729 } 1730 } 1731 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1732 PetscFunctionReturn(0); 1733 } 1734 1735 #undef __FUNCT__ 1736 #define __FUNCT__ "MatGetColumnVector_SeqDense" 1737 PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col) 1738 { 1739 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1740 PetscErrorCode ierr; 1741 PetscScalar *x; 1742 1743 PetscFunctionBegin; 1744 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1745 1746 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1747 ierr = PetscMemcpy(x,a->v+col*a->lda,A->rmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1748 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1749 PetscFunctionReturn(0); 1750 } 1751 1752 /* -------------------------------------------------------------------*/ 1753 static struct _MatOps MatOps_Values = {MatSetValues_SeqDense, 1754 MatGetRow_SeqDense, 1755 MatRestoreRow_SeqDense, 1756 MatMult_SeqDense, 1757 /* 4*/ MatMultAdd_SeqDense, 1758 MatMultTranspose_SeqDense, 1759 MatMultTransposeAdd_SeqDense, 1760 0, 1761 0, 1762 0, 1763 /*10*/ 0, 1764 MatLUFactor_SeqDense, 1765 MatCholeskyFactor_SeqDense, 1766 MatRelax_SeqDense, 1767 MatTranspose_SeqDense, 1768 /*15*/ MatGetInfo_SeqDense, 1769 MatEqual_SeqDense, 1770 MatGetDiagonal_SeqDense, 1771 MatDiagonalScale_SeqDense, 1772 MatNorm_SeqDense, 1773 /*20*/ MatAssemblyBegin_SeqDense, 1774 MatAssemblyEnd_SeqDense, 1775 0, 1776 MatSetOption_SeqDense, 1777 MatZeroEntries_SeqDense, 1778 /*25*/ MatZeroRows_SeqDense, 1779 0, 1780 0, 1781 0, 1782 0, 1783 /*30*/ MatSetUpPreallocation_SeqDense, 1784 0, 1785 0, 1786 MatGetArray_SeqDense, 1787 MatRestoreArray_SeqDense, 1788 /*35*/ MatDuplicate_SeqDense, 1789 0, 1790 0, 1791 0, 1792 0, 1793 /*40*/ MatAXPY_SeqDense, 1794 MatGetSubMatrices_SeqDense, 1795 0, 1796 MatGetValues_SeqDense, 1797 MatCopy_SeqDense, 1798 /*45*/ MatGetRowMax_SeqDense, 1799 MatScale_SeqDense, 1800 0, 1801 0, 1802 0, 1803 /*50*/ 0, 1804 0, 1805 0, 1806 0, 1807 0, 1808 /*55*/ 0, 1809 0, 1810 0, 1811 0, 1812 0, 1813 /*60*/ 0, 1814 MatDestroy_SeqDense, 1815 MatView_SeqDense, 1816 0, 1817 0, 1818 /*65*/ 0, 1819 0, 1820 0, 1821 0, 1822 0, 1823 /*70*/ MatGetRowMaxAbs_SeqDense, 1824 0, 1825 0, 1826 0, 1827 0, 1828 /*75*/ 0, 1829 0, 1830 0, 1831 0, 1832 0, 1833 /*80*/ 0, 1834 0, 1835 0, 1836 0, 1837 /*84*/ MatLoad_SeqDense, 1838 0, 1839 MatIsHermitian_SeqDense, 1840 0, 1841 0, 1842 0, 1843 /*90*/ MatMatMult_SeqDense_SeqDense, 1844 MatMatMultSymbolic_SeqDense_SeqDense, 1845 MatMatMultNumeric_SeqDense_SeqDense, 1846 0, 1847 0, 1848 /*95*/ 0, 1849 MatMatMultTranspose_SeqDense_SeqDense, 1850 MatMatMultTransposeSymbolic_SeqDense_SeqDense, 1851 MatMatMultTransposeNumeric_SeqDense_SeqDense, 1852 0, 1853 /*100*/0, 1854 0, 1855 0, 1856 0, 1857 MatSetSizes_SeqDense, 1858 0, 1859 0, 1860 0, 1861 0, 1862 0, 1863 /*110*/0, 1864 0, 1865 MatGetRowMin_SeqDense, 1866 MatGetColumnVector_SeqDense 1867 }; 1868 1869 #undef __FUNCT__ 1870 #define __FUNCT__ "MatCreateSeqDense" 1871 /*@C 1872 MatCreateSeqDense - Creates a sequential dense matrix that 1873 is stored in column major order (the usual Fortran 77 manner). Many 1874 of the matrix operations use the BLAS and LAPACK routines. 1875 1876 Collective on MPI_Comm 1877 1878 Input Parameters: 1879 + comm - MPI communicator, set to PETSC_COMM_SELF 1880 . m - number of rows 1881 . n - number of columns 1882 - data - optional location of matrix data. Set data=PETSC_NULL for PETSc 1883 to control all matrix memory allocation. 1884 1885 Output Parameter: 1886 . A - the matrix 1887 1888 Notes: 1889 The data input variable is intended primarily for Fortran programmers 1890 who wish to allocate their own matrix memory space. Most users should 1891 set data=PETSC_NULL. 1892 1893 Level: intermediate 1894 1895 .keywords: dense, matrix, LAPACK, BLAS 1896 1897 .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues() 1898 @*/ 1899 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A) 1900 { 1901 PetscErrorCode ierr; 1902 1903 PetscFunctionBegin; 1904 ierr = MatCreate(comm,A);CHKERRQ(ierr); 1905 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 1906 ierr = MatSetType(*A,MATSEQDENSE);CHKERRQ(ierr); 1907 ierr = MatSeqDenseSetPreallocation(*A,data);CHKERRQ(ierr); 1908 PetscFunctionReturn(0); 1909 } 1910 1911 #undef __FUNCT__ 1912 #define __FUNCT__ "MatSeqDenseSetPreallocation" 1913 /*@C 1914 MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements 1915 1916 Collective on MPI_Comm 1917 1918 Input Parameters: 1919 + A - the matrix 1920 - data - the array (or PETSC_NULL) 1921 1922 Notes: 1923 The data input variable is intended primarily for Fortran programmers 1924 who wish to allocate their own matrix memory space. Most users should 1925 need not call this routine. 1926 1927 Level: intermediate 1928 1929 .keywords: dense, matrix, LAPACK, BLAS 1930 1931 .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues() 1932 @*/ 1933 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqDenseSetPreallocation(Mat B,PetscScalar data[]) 1934 { 1935 PetscErrorCode ierr,(*f)(Mat,PetscScalar[]); 1936 1937 PetscFunctionBegin; 1938 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 1939 if (f) { 1940 ierr = (*f)(B,data);CHKERRQ(ierr); 1941 } 1942 PetscFunctionReturn(0); 1943 } 1944 1945 EXTERN_C_BEGIN 1946 #undef __FUNCT__ 1947 #define __FUNCT__ "MatSeqDenseSetPreallocation_SeqDense" 1948 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data) 1949 { 1950 Mat_SeqDense *b; 1951 PetscErrorCode ierr; 1952 1953 PetscFunctionBegin; 1954 B->preallocated = PETSC_TRUE; 1955 b = (Mat_SeqDense*)B->data; 1956 if (b->lda <= 0) b->lda = B->rmap->n; 1957 if (!data) { /* petsc-allocated storage */ 1958 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 1959 ierr = PetscMalloc(b->lda*b->Nmax*sizeof(PetscScalar),&b->v);CHKERRQ(ierr); 1960 ierr = PetscMemzero(b->v,b->lda*b->Nmax*sizeof(PetscScalar));CHKERRQ(ierr); 1961 ierr = PetscLogObjectMemory(B,b->lda*b->Nmax*sizeof(PetscScalar));CHKERRQ(ierr); 1962 b->user_alloc = PETSC_FALSE; 1963 } else { /* user-allocated storage */ 1964 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 1965 b->v = data; 1966 b->user_alloc = PETSC_TRUE; 1967 } 1968 PetscFunctionReturn(0); 1969 } 1970 EXTERN_C_END 1971 1972 #undef __FUNCT__ 1973 #define __FUNCT__ "MatSeqDenseSetLDA" 1974 /*@C 1975 MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array 1976 1977 Input parameter: 1978 + A - the matrix 1979 - lda - the leading dimension 1980 1981 Notes: 1982 This routine is to be used in conjunction with MatSeqDenseSetPreallocation; 1983 it asserts that the preallocation has a leading dimension (the LDA parameter 1984 of Blas and Lapack fame) larger than M, the first dimension of the matrix. 1985 1986 Level: intermediate 1987 1988 .keywords: dense, matrix, LAPACK, BLAS 1989 1990 .seealso: MatCreate(), MatCreateSeqDense(), MatSeqDenseSetPreallocation(), MatSetMaximumSize() 1991 @*/ 1992 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqDenseSetLDA(Mat B,PetscInt lda) 1993 { 1994 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1995 1996 PetscFunctionBegin; 1997 if (lda < B->rmap->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"LDA %D must be at least matrix dimension %D",lda,B->rmap->n); 1998 b->lda = lda; 1999 b->changelda = PETSC_FALSE; 2000 b->Mmax = PetscMax(b->Mmax,lda); 2001 PetscFunctionReturn(0); 2002 } 2003 2004 /*MC 2005 MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices. 2006 2007 Options Database Keys: 2008 . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions() 2009 2010 Level: beginner 2011 2012 .seealso: MatCreateSeqDense() 2013 2014 M*/ 2015 2016 EXTERN_C_BEGIN 2017 #undef __FUNCT__ 2018 #define __FUNCT__ "MatCreate_SeqDense" 2019 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_SeqDense(Mat B) 2020 { 2021 Mat_SeqDense *b; 2022 PetscErrorCode ierr; 2023 PetscMPIInt size; 2024 2025 PetscFunctionBegin; 2026 ierr = MPI_Comm_size(((PetscObject)B)->comm,&size);CHKERRQ(ierr); 2027 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 2028 2029 B->rmap->bs = B->cmap->bs = 1; 2030 ierr = PetscMapSetUp(B->rmap);CHKERRQ(ierr); 2031 ierr = PetscMapSetUp(B->cmap);CHKERRQ(ierr); 2032 2033 ierr = PetscNewLog(B,Mat_SeqDense,&b);CHKERRQ(ierr); 2034 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2035 B->mapping = 0; 2036 B->data = (void*)b; 2037 2038 2039 b->pivots = 0; 2040 b->roworiented = PETSC_TRUE; 2041 b->v = 0; 2042 b->lda = B->rmap->n; 2043 b->changelda = PETSC_FALSE; 2044 b->Mmax = B->rmap->n; 2045 b->Nmax = B->cmap->n; 2046 2047 2048 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_seqdense_petsc_C", 2049 "MatGetFactor_seqdense_petsc", 2050 MatGetFactor_seqdense_petsc);CHKERRQ(ierr); 2051 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqDenseSetPreallocation_C", 2052 "MatSeqDenseSetPreallocation_SeqDense", 2053 MatSeqDenseSetPreallocation_SeqDense);CHKERRQ(ierr); 2054 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMult_seqaij_seqdense_C", 2055 "MatMatMult_SeqAIJ_SeqDense", 2056 MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2057 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C", 2058 "MatMatMultSymbolic_SeqAIJ_SeqDense", 2059 MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2060 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C", 2061 "MatMatMultNumeric_SeqAIJ_SeqDense", 2062 MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2063 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);CHKERRQ(ierr); 2064 PetscFunctionReturn(0); 2065 } 2066 2067 2068 EXTERN_C_END 2069