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