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