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