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