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