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