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 = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 820 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 821 if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 822 PetscFunctionReturn(0); /* do nothing for now */ 823 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 824 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 825 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)A,viewer,"Matrix Object");CHKERRQ(ierr); 826 for (i=0; i<A->rmap->n; i++) { 827 v = a->v + i; 828 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 829 for (j=0; j<A->cmap->n; j++) { 830 #if defined(PETSC_USE_COMPLEX) 831 if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) { 832 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i) ",j,PetscRealPart(*v),PetscImaginaryPart(*v));CHKERRQ(ierr); 833 } else if (PetscRealPart(*v)) { 834 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",j,PetscRealPart(*v));CHKERRQ(ierr); 835 } 836 #else 837 if (*v) { 838 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",j,*v);CHKERRQ(ierr); 839 } 840 #endif 841 v += a->lda; 842 } 843 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 844 } 845 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 846 } else { 847 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 848 #if defined(PETSC_USE_COMPLEX) 849 /* determine if matrix has all real values */ 850 v = a->v; 851 for (i=0; i<A->rmap->n*A->cmap->n; i++) { 852 if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break ;} 853 } 854 #endif 855 if (format == PETSC_VIEWER_ASCII_MATLAB) { 856 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 857 ierr = PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",A->rmap->n,A->cmap->n);CHKERRQ(ierr); 858 ierr = PetscViewerASCIIPrintf(viewer,"%s = zeros(%D,%D);\n",name,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 859 ierr = PetscViewerASCIIPrintf(viewer,"%s = [\n",name);CHKERRQ(ierr); 860 } else { 861 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)A,viewer,"Matrix Object");CHKERRQ(ierr); 862 } 863 864 for (i=0; i<A->rmap->n; i++) { 865 v = a->v + i; 866 for (j=0; j<A->cmap->n; j++) { 867 #if defined(PETSC_USE_COMPLEX) 868 if (allreal) { 869 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e ",PetscRealPart(*v));CHKERRQ(ierr); 870 } else { 871 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16e i ",PetscRealPart(*v),PetscImaginaryPart(*v));CHKERRQ(ierr); 872 } 873 #else 874 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e ",*v);CHKERRQ(ierr); 875 #endif 876 v += a->lda; 877 } 878 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 879 } 880 if (format == PETSC_VIEWER_ASCII_MATLAB) { 881 ierr = PetscViewerASCIIPrintf(viewer,"];\n");CHKERRQ(ierr); 882 } 883 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 884 } 885 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 886 PetscFunctionReturn(0); 887 } 888 889 #undef __FUNCT__ 890 #define __FUNCT__ "MatView_SeqDense_Binary" 891 static PetscErrorCode MatView_SeqDense_Binary(Mat A,PetscViewer viewer) 892 { 893 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 894 PetscErrorCode ierr; 895 int fd; 896 PetscInt ict,j,n = A->cmap->n,m = A->rmap->n,i,*col_lens,nz = m*n; 897 PetscScalar *v,*anonz,*vals; 898 PetscViewerFormat format; 899 900 PetscFunctionBegin; 901 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 902 903 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 904 if (format == PETSC_VIEWER_NATIVE) { 905 /* store the matrix as a dense matrix */ 906 ierr = PetscMalloc(4*sizeof(PetscInt),&col_lens);CHKERRQ(ierr); 907 col_lens[0] = MAT_FILE_CLASSID; 908 col_lens[1] = m; 909 col_lens[2] = n; 910 col_lens[3] = MATRIX_BINARY_FORMAT_DENSE; 911 ierr = PetscBinaryWrite(fd,col_lens,4,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 912 ierr = PetscFree(col_lens);CHKERRQ(ierr); 913 914 /* write out matrix, by rows */ 915 ierr = PetscMalloc((m*n+1)*sizeof(PetscScalar),&vals);CHKERRQ(ierr); 916 v = a->v; 917 for (j=0; j<n; j++) { 918 for (i=0; i<m; i++) { 919 vals[j + i*n] = *v++; 920 } 921 } 922 ierr = PetscBinaryWrite(fd,vals,n*m,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 923 ierr = PetscFree(vals);CHKERRQ(ierr); 924 } else { 925 ierr = PetscMalloc((4+nz)*sizeof(PetscInt),&col_lens);CHKERRQ(ierr); 926 col_lens[0] = MAT_FILE_CLASSID; 927 col_lens[1] = m; 928 col_lens[2] = n; 929 col_lens[3] = nz; 930 931 /* store lengths of each row and write (including header) to file */ 932 for (i=0; i<m; i++) col_lens[4+i] = n; 933 ierr = PetscBinaryWrite(fd,col_lens,4+m,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 934 935 /* Possibly should write in smaller increments, not whole matrix at once? */ 936 /* store column indices (zero start index) */ 937 ict = 0; 938 for (i=0; i<m; i++) { 939 for (j=0; j<n; j++) col_lens[ict++] = j; 940 } 941 ierr = PetscBinaryWrite(fd,col_lens,nz,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 942 ierr = PetscFree(col_lens);CHKERRQ(ierr); 943 944 /* store nonzero values */ 945 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&anonz);CHKERRQ(ierr); 946 ict = 0; 947 for (i=0; i<m; i++) { 948 v = a->v + i; 949 for (j=0; j<n; j++) { 950 anonz[ict++] = *v; v += a->lda; 951 } 952 } 953 ierr = PetscBinaryWrite(fd,anonz,nz,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 954 ierr = PetscFree(anonz);CHKERRQ(ierr); 955 } 956 PetscFunctionReturn(0); 957 } 958 959 #undef __FUNCT__ 960 #define __FUNCT__ "MatView_SeqDense_Draw_Zoom" 961 PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa) 962 { 963 Mat A = (Mat) Aa; 964 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 965 PetscErrorCode ierr; 966 PetscInt m = A->rmap->n,n = A->cmap->n,color,i,j; 967 PetscScalar *v = a->v; 968 PetscViewer viewer; 969 PetscDraw popup; 970 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r,scale,maxv = 0.0; 971 PetscViewerFormat format; 972 973 PetscFunctionBegin; 974 975 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 976 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 977 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 978 979 /* Loop over matrix elements drawing boxes */ 980 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 981 /* Blue for negative and Red for positive */ 982 color = PETSC_DRAW_BLUE; 983 for(j = 0; j < n; j++) { 984 x_l = j; 985 x_r = x_l + 1.0; 986 for(i = 0; i < m; i++) { 987 y_l = m - i - 1.0; 988 y_r = y_l + 1.0; 989 #if defined(PETSC_USE_COMPLEX) 990 if (PetscRealPart(v[j*m+i]) > 0.) { 991 color = PETSC_DRAW_RED; 992 } else if (PetscRealPart(v[j*m+i]) < 0.) { 993 color = PETSC_DRAW_BLUE; 994 } else { 995 continue; 996 } 997 #else 998 if (v[j*m+i] > 0.) { 999 color = PETSC_DRAW_RED; 1000 } else if (v[j*m+i] < 0.) { 1001 color = PETSC_DRAW_BLUE; 1002 } else { 1003 continue; 1004 } 1005 #endif 1006 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 1007 } 1008 } 1009 } else { 1010 /* use contour shading to indicate magnitude of values */ 1011 /* first determine max of all nonzero values */ 1012 for(i = 0; i < m*n; i++) { 1013 if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]); 1014 } 1015 scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/maxv; 1016 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 1017 if (popup) {ierr = PetscDrawScalePopup(popup,0.0,maxv);CHKERRQ(ierr);} 1018 for(j = 0; j < n; j++) { 1019 x_l = j; 1020 x_r = x_l + 1.0; 1021 for(i = 0; i < m; i++) { 1022 y_l = m - i - 1.0; 1023 y_r = y_l + 1.0; 1024 color = PETSC_DRAW_BASIC_COLORS + (int)(scale*PetscAbsScalar(v[j*m+i])); 1025 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 1026 } 1027 } 1028 } 1029 PetscFunctionReturn(0); 1030 } 1031 1032 #undef __FUNCT__ 1033 #define __FUNCT__ "MatView_SeqDense_Draw" 1034 PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer) 1035 { 1036 PetscDraw draw; 1037 PetscBool isnull; 1038 PetscReal xr,yr,xl,yl,h,w; 1039 PetscErrorCode ierr; 1040 1041 PetscFunctionBegin; 1042 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 1043 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 1044 if (isnull) PetscFunctionReturn(0); 1045 1046 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 1047 xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0; 1048 xr += w; yr += h; xl = -w; yl = -h; 1049 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 1050 ierr = PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);CHKERRQ(ierr); 1051 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);CHKERRQ(ierr); 1052 PetscFunctionReturn(0); 1053 } 1054 1055 #undef __FUNCT__ 1056 #define __FUNCT__ "MatView_SeqDense" 1057 PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer) 1058 { 1059 PetscErrorCode ierr; 1060 PetscBool iascii,isbinary,isdraw; 1061 1062 PetscFunctionBegin; 1063 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 1064 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);CHKERRQ(ierr); 1065 ierr = PetscTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);CHKERRQ(ierr); 1066 1067 if (iascii) { 1068 ierr = MatView_SeqDense_ASCII(A,viewer);CHKERRQ(ierr); 1069 } else if (isbinary) { 1070 ierr = MatView_SeqDense_Binary(A,viewer);CHKERRQ(ierr); 1071 } else if (isdraw) { 1072 ierr = MatView_SeqDense_Draw(A,viewer);CHKERRQ(ierr); 1073 } else { 1074 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Viewer type %s not supported by dense matrix",((PetscObject)viewer)->type_name); 1075 } 1076 PetscFunctionReturn(0); 1077 } 1078 1079 #undef __FUNCT__ 1080 #define __FUNCT__ "MatDestroy_SeqDense" 1081 PetscErrorCode MatDestroy_SeqDense(Mat mat) 1082 { 1083 Mat_SeqDense *l = (Mat_SeqDense*)mat->data; 1084 PetscErrorCode ierr; 1085 1086 PetscFunctionBegin; 1087 #if defined(PETSC_USE_LOG) 1088 PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->rmap->n,mat->cmap->n); 1089 #endif 1090 ierr = PetscFree(l->pivots);CHKERRQ(ierr); 1091 if (!l->user_alloc) {ierr = PetscFree(l->v);CHKERRQ(ierr);} 1092 ierr = PetscFree(l);CHKERRQ(ierr); 1093 1094 ierr = PetscObjectChangeTypeName((PetscObject)mat,0);CHKERRQ(ierr); 1095 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatSeqDenseSetPreallocation_C","",PETSC_NULL);CHKERRQ(ierr); 1096 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatMatMult_seqaij_seqdense_C","",PETSC_NULL);CHKERRQ(ierr); 1097 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatMatMultSymbolic_seqaij_seqdense_C","",PETSC_NULL);CHKERRQ(ierr); 1098 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatMatMultNumeric_seqaij_seqdense_C","",PETSC_NULL);CHKERRQ(ierr); 1099 PetscFunctionReturn(0); 1100 } 1101 1102 #undef __FUNCT__ 1103 #define __FUNCT__ "MatTranspose_SeqDense" 1104 PetscErrorCode MatTranspose_SeqDense(Mat A,MatReuse reuse,Mat *matout) 1105 { 1106 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1107 PetscErrorCode ierr; 1108 PetscInt k,j,m,n,M; 1109 PetscScalar *v,tmp; 1110 1111 PetscFunctionBegin; 1112 v = mat->v; m = A->rmap->n; M = mat->lda; n = A->cmap->n; 1113 if (reuse == MAT_REUSE_MATRIX && *matout == A) { /* in place transpose */ 1114 if (m != n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Can not transpose non-square matrix in place"); 1115 else { 1116 for (j=0; j<m; j++) { 1117 for (k=0; k<j; k++) { 1118 tmp = v[j + k*M]; 1119 v[j + k*M] = v[k + j*M]; 1120 v[k + j*M] = tmp; 1121 } 1122 } 1123 } 1124 } else { /* out-of-place transpose */ 1125 Mat tmat; 1126 Mat_SeqDense *tmatd; 1127 PetscScalar *v2; 1128 1129 if (reuse == MAT_INITIAL_MATRIX) { 1130 ierr = MatCreate(((PetscObject)A)->comm,&tmat);CHKERRQ(ierr); 1131 ierr = MatSetSizes(tmat,A->cmap->n,A->rmap->n,A->cmap->n,A->rmap->n);CHKERRQ(ierr); 1132 ierr = MatSetType(tmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1133 ierr = MatSeqDenseSetPreallocation(tmat,PETSC_NULL);CHKERRQ(ierr); 1134 } else { 1135 tmat = *matout; 1136 } 1137 tmatd = (Mat_SeqDense*)tmat->data; 1138 v = mat->v; v2 = tmatd->v; 1139 for (j=0; j<n; j++) { 1140 for (k=0; k<m; k++) v2[j + k*n] = v[k + j*M]; 1141 } 1142 ierr = MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1143 ierr = MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1144 *matout = tmat; 1145 } 1146 PetscFunctionReturn(0); 1147 } 1148 1149 #undef __FUNCT__ 1150 #define __FUNCT__ "MatEqual_SeqDense" 1151 PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscBool *flg) 1152 { 1153 Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data; 1154 Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data; 1155 PetscInt i,j; 1156 PetscScalar *v1 = mat1->v,*v2 = mat2->v; 1157 1158 PetscFunctionBegin; 1159 if (A1->rmap->n != A2->rmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1160 if (A1->cmap->n != A2->cmap->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1161 for (i=0; i<A1->rmap->n; i++) { 1162 v1 = mat1->v+i; v2 = mat2->v+i; 1163 for (j=0; j<A1->cmap->n; j++) { 1164 if (*v1 != *v2) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1165 v1 += mat1->lda; v2 += mat2->lda; 1166 } 1167 } 1168 *flg = PETSC_TRUE; 1169 PetscFunctionReturn(0); 1170 } 1171 1172 #undef __FUNCT__ 1173 #define __FUNCT__ "MatGetDiagonal_SeqDense" 1174 PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v) 1175 { 1176 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1177 PetscErrorCode ierr; 1178 PetscInt i,n,len; 1179 PetscScalar *x,zero = 0.0; 1180 1181 PetscFunctionBegin; 1182 ierr = VecSet(v,zero);CHKERRQ(ierr); 1183 ierr = VecGetSize(v,&n);CHKERRQ(ierr); 1184 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1185 len = PetscMin(A->rmap->n,A->cmap->n); 1186 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec"); 1187 for (i=0; i<len; i++) { 1188 x[i] = mat->v[i*mat->lda + i]; 1189 } 1190 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1191 PetscFunctionReturn(0); 1192 } 1193 1194 #undef __FUNCT__ 1195 #define __FUNCT__ "MatDiagonalScale_SeqDense" 1196 PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr) 1197 { 1198 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1199 PetscScalar *l,*r,x,*v; 1200 PetscErrorCode ierr; 1201 PetscInt i,j,m = A->rmap->n,n = A->cmap->n; 1202 1203 PetscFunctionBegin; 1204 if (ll) { 1205 ierr = VecGetSize(ll,&m);CHKERRQ(ierr); 1206 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1207 if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size"); 1208 for (i=0; i<m; i++) { 1209 x = l[i]; 1210 v = mat->v + i; 1211 for (j=0; j<n; j++) { (*v) *= x; v+= m;} 1212 } 1213 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1214 ierr = PetscLogFlops(n*m);CHKERRQ(ierr); 1215 } 1216 if (rr) { 1217 ierr = VecGetSize(rr,&n);CHKERRQ(ierr); 1218 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1219 if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size"); 1220 for (i=0; i<n; i++) { 1221 x = r[i]; 1222 v = mat->v + i*m; 1223 for (j=0; j<m; j++) { (*v++) *= x;} 1224 } 1225 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1226 ierr = PetscLogFlops(n*m);CHKERRQ(ierr); 1227 } 1228 PetscFunctionReturn(0); 1229 } 1230 1231 #undef __FUNCT__ 1232 #define __FUNCT__ "MatNorm_SeqDense" 1233 PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm) 1234 { 1235 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1236 PetscScalar *v = mat->v; 1237 PetscReal sum = 0.0; 1238 PetscInt lda=mat->lda,m=A->rmap->n,i,j; 1239 PetscErrorCode ierr; 1240 1241 PetscFunctionBegin; 1242 if (type == NORM_FROBENIUS) { 1243 if (lda>m) { 1244 for (j=0; j<A->cmap->n; j++) { 1245 v = mat->v+j*lda; 1246 for (i=0; i<m; i++) { 1247 #if defined(PETSC_USE_COMPLEX) 1248 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1249 #else 1250 sum += (*v)*(*v); v++; 1251 #endif 1252 } 1253 } 1254 } else { 1255 for (i=0; i<A->cmap->n*A->rmap->n; i++) { 1256 #if defined(PETSC_USE_COMPLEX) 1257 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1258 #else 1259 sum += (*v)*(*v); v++; 1260 #endif 1261 } 1262 } 1263 *nrm = sqrt(sum); 1264 ierr = PetscLogFlops(2.0*A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1265 } else if (type == NORM_1) { 1266 *nrm = 0.0; 1267 for (j=0; j<A->cmap->n; j++) { 1268 v = mat->v + j*mat->lda; 1269 sum = 0.0; 1270 for (i=0; i<A->rmap->n; i++) { 1271 sum += PetscAbsScalar(*v); v++; 1272 } 1273 if (sum > *nrm) *nrm = sum; 1274 } 1275 ierr = PetscLogFlops(A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1276 } else if (type == NORM_INFINITY) { 1277 *nrm = 0.0; 1278 for (j=0; j<A->rmap->n; j++) { 1279 v = mat->v + j; 1280 sum = 0.0; 1281 for (i=0; i<A->cmap->n; i++) { 1282 sum += PetscAbsScalar(*v); v += mat->lda; 1283 } 1284 if (sum > *nrm) *nrm = sum; 1285 } 1286 ierr = PetscLogFlops(A->cmap->n*A->rmap->n);CHKERRQ(ierr); 1287 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm"); 1288 PetscFunctionReturn(0); 1289 } 1290 1291 #undef __FUNCT__ 1292 #define __FUNCT__ "MatSetOption_SeqDense" 1293 PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscBool flg) 1294 { 1295 Mat_SeqDense *aij = (Mat_SeqDense*)A->data; 1296 PetscErrorCode ierr; 1297 1298 PetscFunctionBegin; 1299 switch (op) { 1300 case MAT_ROW_ORIENTED: 1301 aij->roworiented = flg; 1302 break; 1303 case MAT_NEW_NONZERO_LOCATIONS: 1304 case MAT_NEW_NONZERO_LOCATION_ERR: 1305 case MAT_NEW_NONZERO_ALLOCATION_ERR: 1306 case MAT_NEW_DIAGONALS: 1307 case MAT_IGNORE_OFF_PROC_ENTRIES: 1308 case MAT_USE_HASH_TABLE: 1309 case MAT_SYMMETRIC: 1310 case MAT_STRUCTURALLY_SYMMETRIC: 1311 case MAT_HERMITIAN: 1312 case MAT_SYMMETRY_ETERNAL: 1313 case MAT_IGNORE_LOWER_TRIANGULAR: 1314 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 1315 break; 1316 default: 1317 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %s",MatOptions[op]); 1318 } 1319 PetscFunctionReturn(0); 1320 } 1321 1322 #undef __FUNCT__ 1323 #define __FUNCT__ "MatZeroEntries_SeqDense" 1324 PetscErrorCode MatZeroEntries_SeqDense(Mat A) 1325 { 1326 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1327 PetscErrorCode ierr; 1328 PetscInt lda=l->lda,m=A->rmap->n,j; 1329 1330 PetscFunctionBegin; 1331 if (lda>m) { 1332 for (j=0; j<A->cmap->n; j++) { 1333 ierr = PetscMemzero(l->v+j*lda,m*sizeof(PetscScalar));CHKERRQ(ierr); 1334 } 1335 } else { 1336 ierr = PetscMemzero(l->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1337 } 1338 PetscFunctionReturn(0); 1339 } 1340 1341 #undef __FUNCT__ 1342 #define __FUNCT__ "MatZeroRows_SeqDense" 1343 PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag) 1344 { 1345 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1346 PetscInt n = A->cmap->n,i,j; 1347 PetscScalar *slot; 1348 1349 PetscFunctionBegin; 1350 for (i=0; i<N; i++) { 1351 slot = l->v + rows[i]; 1352 for (j=0; j<n; j++) { *slot = 0.0; slot += n;} 1353 } 1354 if (diag != 0.0) { 1355 for (i=0; i<N; i++) { 1356 slot = l->v + (n+1)*rows[i]; 1357 *slot = diag; 1358 } 1359 } 1360 PetscFunctionReturn(0); 1361 } 1362 1363 #undef __FUNCT__ 1364 #define __FUNCT__ "MatGetArray_SeqDense" 1365 PetscErrorCode MatGetArray_SeqDense(Mat A,PetscScalar *array[]) 1366 { 1367 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1368 1369 PetscFunctionBegin; 1370 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"); 1371 *array = mat->v; 1372 PetscFunctionReturn(0); 1373 } 1374 1375 #undef __FUNCT__ 1376 #define __FUNCT__ "MatRestoreArray_SeqDense" 1377 PetscErrorCode MatRestoreArray_SeqDense(Mat A,PetscScalar *array[]) 1378 { 1379 PetscFunctionBegin; 1380 *array = 0; /* user cannot accidently use the array later */ 1381 PetscFunctionReturn(0); 1382 } 1383 1384 #undef __FUNCT__ 1385 #define __FUNCT__ "MatGetSubMatrix_SeqDense" 1386 static PetscErrorCode MatGetSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B) 1387 { 1388 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1389 PetscErrorCode ierr; 1390 PetscInt i,j,nrows,ncols; 1391 const PetscInt *irow,*icol; 1392 PetscScalar *av,*bv,*v = mat->v; 1393 Mat newmat; 1394 1395 PetscFunctionBegin; 1396 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1397 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1398 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1399 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1400 1401 /* Check submatrixcall */ 1402 if (scall == MAT_REUSE_MATRIX) { 1403 PetscInt n_cols,n_rows; 1404 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1405 if (n_rows != nrows || n_cols != ncols) { 1406 /* resize the result result matrix to match number of requested rows/columns */ 1407 ierr = MatSetSizes(*B,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1408 } 1409 newmat = *B; 1410 } else { 1411 /* Create and fill new matrix */ 1412 ierr = MatCreate(((PetscObject)A)->comm,&newmat);CHKERRQ(ierr); 1413 ierr = MatSetSizes(newmat,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1414 ierr = MatSetType(newmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1415 ierr = MatSeqDenseSetPreallocation(newmat,PETSC_NULL);CHKERRQ(ierr); 1416 } 1417 1418 /* Now extract the data pointers and do the copy,column at a time */ 1419 bv = ((Mat_SeqDense*)newmat->data)->v; 1420 1421 for (i=0; i<ncols; i++) { 1422 av = v + mat->lda*icol[i]; 1423 for (j=0; j<nrows; j++) { 1424 *bv++ = av[irow[j]]; 1425 } 1426 } 1427 1428 /* Assemble the matrices so that the correct flags are set */ 1429 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1430 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1431 1432 /* Free work space */ 1433 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1434 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1435 *B = newmat; 1436 PetscFunctionReturn(0); 1437 } 1438 1439 #undef __FUNCT__ 1440 #define __FUNCT__ "MatGetSubMatrices_SeqDense" 1441 PetscErrorCode MatGetSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1442 { 1443 PetscErrorCode ierr; 1444 PetscInt i; 1445 1446 PetscFunctionBegin; 1447 if (scall == MAT_INITIAL_MATRIX) { 1448 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1449 } 1450 1451 for (i=0; i<n; i++) { 1452 ierr = MatGetSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1453 } 1454 PetscFunctionReturn(0); 1455 } 1456 1457 #undef __FUNCT__ 1458 #define __FUNCT__ "MatAssemblyBegin_SeqDense" 1459 PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode) 1460 { 1461 PetscFunctionBegin; 1462 PetscFunctionReturn(0); 1463 } 1464 1465 #undef __FUNCT__ 1466 #define __FUNCT__ "MatAssemblyEnd_SeqDense" 1467 PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode) 1468 { 1469 PetscFunctionBegin; 1470 PetscFunctionReturn(0); 1471 } 1472 1473 #undef __FUNCT__ 1474 #define __FUNCT__ "MatCopy_SeqDense" 1475 PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str) 1476 { 1477 Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense *)B->data; 1478 PetscErrorCode ierr; 1479 PetscInt lda1=a->lda,lda2=b->lda, m=A->rmap->n,n=A->cmap->n, j; 1480 1481 PetscFunctionBegin; 1482 /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */ 1483 if (A->ops->copy != B->ops->copy) { 1484 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1485 PetscFunctionReturn(0); 1486 } 1487 if (m != B->rmap->n || n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"size(B) != size(A)"); 1488 if (lda1>m || lda2>m) { 1489 for (j=0; j<n; j++) { 1490 ierr = PetscMemcpy(b->v+j*lda2,a->v+j*lda1,m*sizeof(PetscScalar));CHKERRQ(ierr); 1491 } 1492 } else { 1493 ierr = PetscMemcpy(b->v,a->v,A->rmap->n*A->cmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1494 } 1495 PetscFunctionReturn(0); 1496 } 1497 1498 #undef __FUNCT__ 1499 #define __FUNCT__ "MatSetUpPreallocation_SeqDense" 1500 PetscErrorCode MatSetUpPreallocation_SeqDense(Mat A) 1501 { 1502 PetscErrorCode ierr; 1503 1504 PetscFunctionBegin; 1505 ierr = MatSeqDenseSetPreallocation(A,0);CHKERRQ(ierr); 1506 PetscFunctionReturn(0); 1507 } 1508 1509 #undef __FUNCT__ 1510 #define __FUNCT__ "MatSetSizes_SeqDense" 1511 PetscErrorCode MatSetSizes_SeqDense(Mat A,PetscInt m,PetscInt n,PetscInt M,PetscInt N) 1512 { 1513 PetscFunctionBegin; 1514 /* this will not be called before lda, Mmax, and Nmax have been set */ 1515 m = PetscMax(m,M); 1516 n = PetscMax(n,N); 1517 1518 /* 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); 1519 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); 1520 */ 1521 A->rmap->n = A->rmap->N = m; 1522 A->cmap->n = A->cmap->N = n; 1523 PetscFunctionReturn(0); 1524 } 1525 1526 #undef __FUNCT__ 1527 #define __FUNCT__ "MatConjugate_SeqDense" 1528 static PetscErrorCode MatConjugate_SeqDense(Mat A) 1529 { 1530 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1531 PetscInt i,nz = A->rmap->n*A->cmap->n; 1532 PetscScalar *aa = a->v; 1533 1534 PetscFunctionBegin; 1535 for (i=0; i<nz; i++) aa[i] = PetscConj(aa[i]); 1536 PetscFunctionReturn(0); 1537 } 1538 1539 #undef __FUNCT__ 1540 #define __FUNCT__ "MatRealPart_SeqDense" 1541 static PetscErrorCode MatRealPart_SeqDense(Mat A) 1542 { 1543 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1544 PetscInt i,nz = A->rmap->n*A->cmap->n; 1545 PetscScalar *aa = a->v; 1546 1547 PetscFunctionBegin; 1548 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 1549 PetscFunctionReturn(0); 1550 } 1551 1552 #undef __FUNCT__ 1553 #define __FUNCT__ "MatImaginaryPart_SeqDense" 1554 static PetscErrorCode MatImaginaryPart_SeqDense(Mat A) 1555 { 1556 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1557 PetscInt i,nz = A->rmap->n*A->cmap->n; 1558 PetscScalar *aa = a->v; 1559 1560 PetscFunctionBegin; 1561 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 1562 PetscFunctionReturn(0); 1563 } 1564 1565 /* ----------------------------------------------------------------*/ 1566 #undef __FUNCT__ 1567 #define __FUNCT__ "MatMatMult_SeqDense_SeqDense" 1568 PetscErrorCode MatMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 1569 { 1570 PetscErrorCode ierr; 1571 1572 PetscFunctionBegin; 1573 if (scall == MAT_INITIAL_MATRIX){ 1574 ierr = MatMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 1575 } 1576 ierr = MatMatMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 1577 PetscFunctionReturn(0); 1578 } 1579 1580 #undef __FUNCT__ 1581 #define __FUNCT__ "MatMatMultSymbolic_SeqDense_SeqDense" 1582 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 1583 { 1584 PetscErrorCode ierr; 1585 PetscInt m=A->rmap->n,n=B->cmap->n; 1586 Mat Cmat; 1587 1588 PetscFunctionBegin; 1589 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); 1590 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 1591 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 1592 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 1593 ierr = MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);CHKERRQ(ierr); 1594 Cmat->assembled = PETSC_TRUE; 1595 *C = Cmat; 1596 PetscFunctionReturn(0); 1597 } 1598 1599 #undef __FUNCT__ 1600 #define __FUNCT__ "MatMatMultNumeric_SeqDense_SeqDense" 1601 PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 1602 { 1603 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1604 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1605 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 1606 PetscBLASInt m,n,k; 1607 PetscScalar _DOne=1.0,_DZero=0.0; 1608 1609 PetscFunctionBegin; 1610 m = PetscBLASIntCast(A->rmap->n); 1611 n = PetscBLASIntCast(B->cmap->n); 1612 k = PetscBLASIntCast(A->cmap->n); 1613 BLASgemm_("N","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda); 1614 PetscFunctionReturn(0); 1615 } 1616 1617 #undef __FUNCT__ 1618 #define __FUNCT__ "MatMatMultTranspose_SeqDense_SeqDense" 1619 PetscErrorCode MatMatMultTranspose_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 1620 { 1621 PetscErrorCode ierr; 1622 1623 PetscFunctionBegin; 1624 if (scall == MAT_INITIAL_MATRIX){ 1625 ierr = MatMatMultTransposeSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 1626 } 1627 ierr = MatMatMultTransposeNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 1628 PetscFunctionReturn(0); 1629 } 1630 1631 #undef __FUNCT__ 1632 #define __FUNCT__ "MatMatMultTransposeSymbolic_SeqDense_SeqDense" 1633 PetscErrorCode MatMatMultTransposeSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 1634 { 1635 PetscErrorCode ierr; 1636 PetscInt m=A->cmap->n,n=B->cmap->n; 1637 Mat Cmat; 1638 1639 PetscFunctionBegin; 1640 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); 1641 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 1642 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 1643 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 1644 ierr = MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);CHKERRQ(ierr); 1645 Cmat->assembled = PETSC_TRUE; 1646 *C = Cmat; 1647 PetscFunctionReturn(0); 1648 } 1649 1650 #undef __FUNCT__ 1651 #define __FUNCT__ "MatMatMultTransposeNumeric_SeqDense_SeqDense" 1652 PetscErrorCode MatMatMultTransposeNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 1653 { 1654 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1655 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1656 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 1657 PetscBLASInt m,n,k; 1658 PetscScalar _DOne=1.0,_DZero=0.0; 1659 1660 PetscFunctionBegin; 1661 m = PetscBLASIntCast(A->cmap->n); 1662 n = PetscBLASIntCast(B->cmap->n); 1663 k = PetscBLASIntCast(A->rmap->n); 1664 /* 1665 Note the m and n arguments below are the number rows and columns of A', not A! 1666 */ 1667 BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda); 1668 PetscFunctionReturn(0); 1669 } 1670 1671 #undef __FUNCT__ 1672 #define __FUNCT__ "MatGetRowMax_SeqDense" 1673 PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[]) 1674 { 1675 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1676 PetscErrorCode ierr; 1677 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1678 PetscScalar *x; 1679 MatScalar *aa = a->v; 1680 1681 PetscFunctionBegin; 1682 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1683 1684 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1685 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1686 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1687 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1688 for (i=0; i<m; i++) { 1689 x[i] = aa[i]; if (idx) idx[i] = 0; 1690 for (j=1; j<n; j++){ 1691 if (PetscRealPart(x[i]) < PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;} 1692 } 1693 } 1694 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1695 PetscFunctionReturn(0); 1696 } 1697 1698 #undef __FUNCT__ 1699 #define __FUNCT__ "MatGetRowMaxAbs_SeqDense" 1700 PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[]) 1701 { 1702 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1703 PetscErrorCode ierr; 1704 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1705 PetscScalar *x; 1706 PetscReal atmp; 1707 MatScalar *aa = a->v; 1708 1709 PetscFunctionBegin; 1710 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1711 1712 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1713 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1714 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1715 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1716 for (i=0; i<m; i++) { 1717 x[i] = PetscAbsScalar(aa[i]); 1718 for (j=1; j<n; j++){ 1719 atmp = PetscAbsScalar(aa[i+m*j]); 1720 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;} 1721 } 1722 } 1723 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1724 PetscFunctionReturn(0); 1725 } 1726 1727 #undef __FUNCT__ 1728 #define __FUNCT__ "MatGetRowMin_SeqDense" 1729 PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[]) 1730 { 1731 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1732 PetscErrorCode ierr; 1733 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,p; 1734 PetscScalar *x; 1735 MatScalar *aa = a->v; 1736 1737 PetscFunctionBegin; 1738 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1739 1740 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1741 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1742 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1743 if (p != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1744 for (i=0; i<m; i++) { 1745 x[i] = aa[i]; if (idx) idx[i] = 0; 1746 for (j=1; j<n; j++){ 1747 if (PetscRealPart(x[i]) > PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;} 1748 } 1749 } 1750 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1751 PetscFunctionReturn(0); 1752 } 1753 1754 #undef __FUNCT__ 1755 #define __FUNCT__ "MatGetColumnVector_SeqDense" 1756 PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col) 1757 { 1758 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1759 PetscErrorCode ierr; 1760 PetscScalar *x; 1761 1762 PetscFunctionBegin; 1763 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1764 1765 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1766 ierr = PetscMemcpy(x,a->v+col*a->lda,A->rmap->n*sizeof(PetscScalar));CHKERRQ(ierr); 1767 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1768 PetscFunctionReturn(0); 1769 } 1770 1771 /* -------------------------------------------------------------------*/ 1772 static struct _MatOps MatOps_Values = {MatSetValues_SeqDense, 1773 MatGetRow_SeqDense, 1774 MatRestoreRow_SeqDense, 1775 MatMult_SeqDense, 1776 /* 4*/ MatMultAdd_SeqDense, 1777 MatMultTranspose_SeqDense, 1778 MatMultTransposeAdd_SeqDense, 1779 0, 1780 0, 1781 0, 1782 /*10*/ 0, 1783 MatLUFactor_SeqDense, 1784 MatCholeskyFactor_SeqDense, 1785 MatSOR_SeqDense, 1786 MatTranspose_SeqDense, 1787 /*15*/ MatGetInfo_SeqDense, 1788 MatEqual_SeqDense, 1789 MatGetDiagonal_SeqDense, 1790 MatDiagonalScale_SeqDense, 1791 MatNorm_SeqDense, 1792 /*20*/ MatAssemblyBegin_SeqDense, 1793 MatAssemblyEnd_SeqDense, 1794 MatSetOption_SeqDense, 1795 MatZeroEntries_SeqDense, 1796 /*24*/ MatZeroRows_SeqDense, 1797 0, 1798 0, 1799 0, 1800 0, 1801 /*29*/ MatSetUpPreallocation_SeqDense, 1802 0, 1803 0, 1804 MatGetArray_SeqDense, 1805 MatRestoreArray_SeqDense, 1806 /*34*/ MatDuplicate_SeqDense, 1807 0, 1808 0, 1809 0, 1810 0, 1811 /*39*/ MatAXPY_SeqDense, 1812 MatGetSubMatrices_SeqDense, 1813 0, 1814 MatGetValues_SeqDense, 1815 MatCopy_SeqDense, 1816 /*44*/ MatGetRowMax_SeqDense, 1817 MatScale_SeqDense, 1818 0, 1819 0, 1820 0, 1821 /*49*/ 0, 1822 0, 1823 0, 1824 0, 1825 0, 1826 /*54*/ 0, 1827 0, 1828 0, 1829 0, 1830 0, 1831 /*59*/ 0, 1832 MatDestroy_SeqDense, 1833 MatView_SeqDense, 1834 0, 1835 0, 1836 /*64*/ 0, 1837 0, 1838 0, 1839 0, 1840 0, 1841 /*69*/ MatGetRowMaxAbs_SeqDense, 1842 0, 1843 0, 1844 0, 1845 0, 1846 /*74*/ 0, 1847 0, 1848 0, 1849 0, 1850 0, 1851 /*79*/ 0, 1852 0, 1853 0, 1854 0, 1855 /*83*/ MatLoad_SeqDense, 1856 0, 1857 MatIsHermitian_SeqDense, 1858 0, 1859 0, 1860 0, 1861 /*89*/ MatMatMult_SeqDense_SeqDense, 1862 MatMatMultSymbolic_SeqDense_SeqDense, 1863 MatMatMultNumeric_SeqDense_SeqDense, 1864 0, 1865 0, 1866 /*94*/ 0, 1867 MatMatMultTranspose_SeqDense_SeqDense, 1868 MatMatMultTransposeSymbolic_SeqDense_SeqDense, 1869 MatMatMultTransposeNumeric_SeqDense_SeqDense, 1870 0, 1871 /*99*/ 0, 1872 0, 1873 0, 1874 MatConjugate_SeqDense, 1875 MatSetSizes_SeqDense, 1876 /*104*/0, 1877 MatRealPart_SeqDense, 1878 MatImaginaryPart_SeqDense, 1879 0, 1880 0, 1881 /*109*/0, 1882 0, 1883 MatGetRowMin_SeqDense, 1884 MatGetColumnVector_SeqDense, 1885 0, 1886 /*114*/0, 1887 0, 1888 0, 1889 0, 1890 0, 1891 /*119*/0, 1892 0, 1893 0, 1894 0 1895 }; 1896 1897 #undef __FUNCT__ 1898 #define __FUNCT__ "MatCreateSeqDense" 1899 /*@C 1900 MatCreateSeqDense - Creates a sequential dense matrix that 1901 is stored in column major order (the usual Fortran 77 manner). Many 1902 of the matrix operations use the BLAS and LAPACK routines. 1903 1904 Collective on MPI_Comm 1905 1906 Input Parameters: 1907 + comm - MPI communicator, set to PETSC_COMM_SELF 1908 . m - number of rows 1909 . n - number of columns 1910 - data - optional location of matrix data in column major order. Set data=PETSC_NULL for PETSc 1911 to control all matrix memory allocation. 1912 1913 Output Parameter: 1914 . A - the matrix 1915 1916 Notes: 1917 The data input variable is intended primarily for Fortran programmers 1918 who wish to allocate their own matrix memory space. Most users should 1919 set data=PETSC_NULL. 1920 1921 Level: intermediate 1922 1923 .keywords: dense, matrix, LAPACK, BLAS 1924 1925 .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues() 1926 @*/ 1927 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A) 1928 { 1929 PetscErrorCode ierr; 1930 1931 PetscFunctionBegin; 1932 ierr = MatCreate(comm,A);CHKERRQ(ierr); 1933 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 1934 ierr = MatSetType(*A,MATSEQDENSE);CHKERRQ(ierr); 1935 ierr = MatSeqDenseSetPreallocation(*A,data);CHKERRQ(ierr); 1936 PetscFunctionReturn(0); 1937 } 1938 1939 #undef __FUNCT__ 1940 #define __FUNCT__ "MatSeqDenseSetPreallocation" 1941 /*@C 1942 MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements 1943 1944 Collective on MPI_Comm 1945 1946 Input Parameters: 1947 + A - the matrix 1948 - data - the array (or PETSC_NULL) 1949 1950 Notes: 1951 The data input variable is intended primarily for Fortran programmers 1952 who wish to allocate their own matrix memory space. Most users should 1953 need not call this routine. 1954 1955 Level: intermediate 1956 1957 .keywords: dense, matrix, LAPACK, BLAS 1958 1959 .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues(), MatSeqDenseSetLDA() 1960 1961 @*/ 1962 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqDenseSetPreallocation(Mat B,PetscScalar data[]) 1963 { 1964 PetscErrorCode ierr,(*f)(Mat,PetscScalar[]); 1965 1966 PetscFunctionBegin; 1967 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 1968 if (f) { 1969 ierr = (*f)(B,data);CHKERRQ(ierr); 1970 } 1971 PetscFunctionReturn(0); 1972 } 1973 1974 EXTERN_C_BEGIN 1975 #undef __FUNCT__ 1976 #define __FUNCT__ "MatSeqDenseSetPreallocation_SeqDense" 1977 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data) 1978 { 1979 Mat_SeqDense *b; 1980 PetscErrorCode ierr; 1981 1982 PetscFunctionBegin; 1983 B->preallocated = PETSC_TRUE; 1984 1985 ierr = PetscLayoutSetBlockSize(B->rmap,1);CHKERRQ(ierr); 1986 ierr = PetscLayoutSetBlockSize(B->cmap,1);CHKERRQ(ierr); 1987 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 1988 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 1989 1990 b = (Mat_SeqDense*)B->data; 1991 b->Mmax = B->rmap->n; 1992 b->Nmax = B->cmap->n; 1993 if(b->lda <= 0 || b->changelda) b->lda = B->rmap->n; 1994 1995 if (!data) { /* petsc-allocated storage */ 1996 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 1997 ierr = PetscMalloc(b->lda*b->Nmax*sizeof(PetscScalar),&b->v);CHKERRQ(ierr); 1998 ierr = PetscMemzero(b->v,b->lda*b->Nmax*sizeof(PetscScalar));CHKERRQ(ierr); 1999 ierr = PetscLogObjectMemory(B,b->lda*b->Nmax*sizeof(PetscScalar));CHKERRQ(ierr); 2000 b->user_alloc = PETSC_FALSE; 2001 } else { /* user-allocated storage */ 2002 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 2003 b->v = data; 2004 b->user_alloc = PETSC_TRUE; 2005 } 2006 B->assembled = PETSC_TRUE; 2007 PetscFunctionReturn(0); 2008 } 2009 EXTERN_C_END 2010 2011 #undef __FUNCT__ 2012 #define __FUNCT__ "MatSeqDenseSetLDA" 2013 /*@C 2014 MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array 2015 2016 Input parameter: 2017 + A - the matrix 2018 - lda - the leading dimension 2019 2020 Notes: 2021 This routine is to be used in conjunction with MatSeqDenseSetPreallocation(); 2022 it asserts that the preallocation has a leading dimension (the LDA parameter 2023 of Blas and Lapack fame) larger than M, the first dimension of the matrix. 2024 2025 Level: intermediate 2026 2027 .keywords: dense, matrix, LAPACK, BLAS 2028 2029 .seealso: MatCreate(), MatCreateSeqDense(), MatSeqDenseSetPreallocation(), MatSetMaximumSize() 2030 2031 @*/ 2032 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqDenseSetLDA(Mat B,PetscInt lda) 2033 { 2034 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 2035 2036 PetscFunctionBegin; 2037 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); 2038 b->lda = lda; 2039 b->changelda = PETSC_FALSE; 2040 b->Mmax = PetscMax(b->Mmax,lda); 2041 PetscFunctionReturn(0); 2042 } 2043 2044 /*MC 2045 MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices. 2046 2047 Options Database Keys: 2048 . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions() 2049 2050 Level: beginner 2051 2052 .seealso: MatCreateSeqDense() 2053 2054 M*/ 2055 2056 EXTERN_C_BEGIN 2057 #undef __FUNCT__ 2058 #define __FUNCT__ "MatCreate_SeqDense" 2059 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_SeqDense(Mat B) 2060 { 2061 Mat_SeqDense *b; 2062 PetscErrorCode ierr; 2063 PetscMPIInt size; 2064 2065 PetscFunctionBegin; 2066 ierr = MPI_Comm_size(((PetscObject)B)->comm,&size);CHKERRQ(ierr); 2067 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 2068 2069 ierr = PetscNewLog(B,Mat_SeqDense,&b);CHKERRQ(ierr); 2070 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2071 B->mapping = 0; 2072 B->data = (void*)b; 2073 2074 b->pivots = 0; 2075 b->roworiented = PETSC_TRUE; 2076 b->v = 0; 2077 b->changelda = PETSC_FALSE; 2078 2079 2080 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_petsc_C", 2081 "MatGetFactor_seqdense_petsc", 2082 MatGetFactor_seqdense_petsc);CHKERRQ(ierr); 2083 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqDenseSetPreallocation_C", 2084 "MatSeqDenseSetPreallocation_SeqDense", 2085 MatSeqDenseSetPreallocation_SeqDense);CHKERRQ(ierr); 2086 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMult_seqaij_seqdense_C", 2087 "MatMatMult_SeqAIJ_SeqDense", 2088 MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2089 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C", 2090 "MatMatMultSymbolic_SeqAIJ_SeqDense", 2091 MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2092 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C", 2093 "MatMatMultNumeric_SeqAIJ_SeqDense", 2094 MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2095 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);CHKERRQ(ierr); 2096 PetscFunctionReturn(0); 2097 } 2098 EXTERN_C_END 2099