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