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