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