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