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