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