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