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,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) { /* 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 ierr = MatCreate(((PetscObject)A)->comm,&tmat);CHKERRQ(ierr); 1106 ierr = MatSetSizes(tmat,A->cmap.n,A->rmap.n,A->cmap.n,A->rmap.n);CHKERRQ(ierr); 1107 ierr = MatSetType(tmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1108 ierr = MatSeqDenseSetPreallocation(tmat,PETSC_NULL);CHKERRQ(ierr); 1109 tmatd = (Mat_SeqDense*)tmat->data; 1110 v = mat->v; v2 = tmatd->v; 1111 for (j=0; j<n; j++) { 1112 for (k=0; k<m; k++) v2[j + k*n] = v[k + j*M]; 1113 } 1114 ierr = MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1115 ierr = MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1116 *matout = tmat; 1117 } 1118 PetscFunctionReturn(0); 1119 } 1120 1121 #undef __FUNCT__ 1122 #define __FUNCT__ "MatEqual_SeqDense" 1123 PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscTruth *flg) 1124 { 1125 Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data; 1126 Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data; 1127 PetscInt i,j; 1128 PetscScalar *v1 = mat1->v,*v2 = mat2->v; 1129 1130 PetscFunctionBegin; 1131 if (A1->rmap.n != A2->rmap.n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1132 if (A1->cmap.n != A2->cmap.n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1133 for (i=0; i<A1->rmap.n; i++) { 1134 v1 = mat1->v+i; v2 = mat2->v+i; 1135 for (j=0; j<A1->cmap.n; j++) { 1136 if (*v1 != *v2) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1137 v1 += mat1->lda; v2 += mat2->lda; 1138 } 1139 } 1140 *flg = PETSC_TRUE; 1141 PetscFunctionReturn(0); 1142 } 1143 1144 #undef __FUNCT__ 1145 #define __FUNCT__ "MatGetDiagonal_SeqDense" 1146 PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v) 1147 { 1148 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1149 PetscErrorCode ierr; 1150 PetscInt i,n,len; 1151 PetscScalar *x,zero = 0.0; 1152 1153 PetscFunctionBegin; 1154 ierr = VecSet(v,zero);CHKERRQ(ierr); 1155 ierr = VecGetSize(v,&n);CHKERRQ(ierr); 1156 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1157 len = PetscMin(A->rmap.n,A->cmap.n); 1158 if (n != A->rmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec"); 1159 for (i=0; i<len; i++) { 1160 x[i] = mat->v[i*mat->lda + i]; 1161 } 1162 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1163 PetscFunctionReturn(0); 1164 } 1165 1166 #undef __FUNCT__ 1167 #define __FUNCT__ "MatDiagonalScale_SeqDense" 1168 PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr) 1169 { 1170 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1171 PetscScalar *l,*r,x,*v; 1172 PetscErrorCode ierr; 1173 PetscInt i,j,m = A->rmap.n,n = A->cmap.n; 1174 1175 PetscFunctionBegin; 1176 if (ll) { 1177 ierr = VecGetSize(ll,&m);CHKERRQ(ierr); 1178 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1179 if (m != A->rmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size"); 1180 for (i=0; i<m; i++) { 1181 x = l[i]; 1182 v = mat->v + i; 1183 for (j=0; j<n; j++) { (*v) *= x; v+= m;} 1184 } 1185 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1186 ierr = PetscLogFlops(n*m);CHKERRQ(ierr); 1187 } 1188 if (rr) { 1189 ierr = VecGetSize(rr,&n);CHKERRQ(ierr); 1190 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1191 if (n != A->cmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size"); 1192 for (i=0; i<n; i++) { 1193 x = r[i]; 1194 v = mat->v + i*m; 1195 for (j=0; j<m; j++) { (*v++) *= x;} 1196 } 1197 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1198 ierr = PetscLogFlops(n*m);CHKERRQ(ierr); 1199 } 1200 PetscFunctionReturn(0); 1201 } 1202 1203 #undef __FUNCT__ 1204 #define __FUNCT__ "MatNorm_SeqDense" 1205 PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm) 1206 { 1207 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1208 PetscScalar *v = mat->v; 1209 PetscReal sum = 0.0; 1210 PetscInt lda=mat->lda,m=A->rmap.n,i,j; 1211 PetscErrorCode ierr; 1212 1213 PetscFunctionBegin; 1214 if (type == NORM_FROBENIUS) { 1215 if (lda>m) { 1216 for (j=0; j<A->cmap.n; j++) { 1217 v = mat->v+j*lda; 1218 for (i=0; i<m; i++) { 1219 #if defined(PETSC_USE_COMPLEX) 1220 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1221 #else 1222 sum += (*v)*(*v); v++; 1223 #endif 1224 } 1225 } 1226 } else { 1227 for (i=0; i<A->cmap.n*A->rmap.n; i++) { 1228 #if defined(PETSC_USE_COMPLEX) 1229 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1230 #else 1231 sum += (*v)*(*v); v++; 1232 #endif 1233 } 1234 } 1235 *nrm = sqrt(sum); 1236 ierr = PetscLogFlops(2*A->cmap.n*A->rmap.n);CHKERRQ(ierr); 1237 } else if (type == NORM_1) { 1238 *nrm = 0.0; 1239 for (j=0; j<A->cmap.n; j++) { 1240 v = mat->v + j*mat->lda; 1241 sum = 0.0; 1242 for (i=0; i<A->rmap.n; i++) { 1243 sum += PetscAbsScalar(*v); v++; 1244 } 1245 if (sum > *nrm) *nrm = sum; 1246 } 1247 ierr = PetscLogFlops(A->cmap.n*A->rmap.n);CHKERRQ(ierr); 1248 } else if (type == NORM_INFINITY) { 1249 *nrm = 0.0; 1250 for (j=0; j<A->rmap.n; j++) { 1251 v = mat->v + j; 1252 sum = 0.0; 1253 for (i=0; i<A->cmap.n; i++) { 1254 sum += PetscAbsScalar(*v); v += mat->lda; 1255 } 1256 if (sum > *nrm) *nrm = sum; 1257 } 1258 ierr = PetscLogFlops(A->cmap.n*A->rmap.n);CHKERRQ(ierr); 1259 } else { 1260 SETERRQ(PETSC_ERR_SUP,"No two norm"); 1261 } 1262 PetscFunctionReturn(0); 1263 } 1264 1265 #undef __FUNCT__ 1266 #define __FUNCT__ "MatSetOption_SeqDense" 1267 PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscTruth flg) 1268 { 1269 Mat_SeqDense *aij = (Mat_SeqDense*)A->data; 1270 PetscErrorCode ierr; 1271 1272 PetscFunctionBegin; 1273 switch (op) { 1274 case MAT_ROW_ORIENTED: 1275 aij->roworiented = flg; 1276 break; 1277 case MAT_NEW_NONZERO_LOCATIONS: 1278 case MAT_NEW_NONZERO_LOCATION_ERR: 1279 case MAT_NEW_NONZERO_ALLOCATION_ERR: 1280 case MAT_NEW_DIAGONALS: 1281 case MAT_IGNORE_OFF_PROC_ENTRIES: 1282 case MAT_USE_HASH_TABLE: 1283 case MAT_SYMMETRIC: 1284 case MAT_STRUCTURALLY_SYMMETRIC: 1285 case MAT_HERMITIAN: 1286 case MAT_SYMMETRY_ETERNAL: 1287 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 1288 break; 1289 default: 1290 SETERRQ1(PETSC_ERR_SUP,"unknown option %d",op); 1291 } 1292 PetscFunctionReturn(0); 1293 } 1294 1295 #undef __FUNCT__ 1296 #define __FUNCT__ "MatZeroEntries_SeqDense" 1297 PetscErrorCode MatZeroEntries_SeqDense(Mat A) 1298 { 1299 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1300 PetscErrorCode ierr; 1301 PetscInt lda=l->lda,m=A->rmap.n,j; 1302 1303 PetscFunctionBegin; 1304 if (lda>m) { 1305 for (j=0; j<A->cmap.n; j++) { 1306 ierr = PetscMemzero(l->v+j*lda,m*sizeof(PetscScalar));CHKERRQ(ierr); 1307 } 1308 } else { 1309 ierr = PetscMemzero(l->v,A->rmap.n*A->cmap.n*sizeof(PetscScalar));CHKERRQ(ierr); 1310 } 1311 PetscFunctionReturn(0); 1312 } 1313 1314 #undef __FUNCT__ 1315 #define __FUNCT__ "MatZeroRows_SeqDense" 1316 PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag) 1317 { 1318 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1319 PetscInt n = A->cmap.n,i,j; 1320 PetscScalar *slot; 1321 1322 PetscFunctionBegin; 1323 for (i=0; i<N; i++) { 1324 slot = l->v + rows[i]; 1325 for (j=0; j<n; j++) { *slot = 0.0; slot += n;} 1326 } 1327 if (diag != 0.0) { 1328 for (i=0; i<N; i++) { 1329 slot = l->v + (n+1)*rows[i]; 1330 *slot = diag; 1331 } 1332 } 1333 PetscFunctionReturn(0); 1334 } 1335 1336 #undef __FUNCT__ 1337 #define __FUNCT__ "MatGetArray_SeqDense" 1338 PetscErrorCode MatGetArray_SeqDense(Mat A,PetscScalar *array[]) 1339 { 1340 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1341 1342 PetscFunctionBegin; 1343 if (mat->lda != A->rmap.n) SETERRQ(PETSC_ERR_SUP,"Cannot get array for Dense matrices with LDA different from number of rows"); 1344 *array = mat->v; 1345 PetscFunctionReturn(0); 1346 } 1347 1348 #undef __FUNCT__ 1349 #define __FUNCT__ "MatRestoreArray_SeqDense" 1350 PetscErrorCode MatRestoreArray_SeqDense(Mat A,PetscScalar *array[]) 1351 { 1352 PetscFunctionBegin; 1353 *array = 0; /* user cannot accidently use the array later */ 1354 PetscFunctionReturn(0); 1355 } 1356 1357 #undef __FUNCT__ 1358 #define __FUNCT__ "MatGetSubMatrix_SeqDense" 1359 static PetscErrorCode MatGetSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B) 1360 { 1361 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1362 PetscErrorCode ierr; 1363 PetscInt i,j,*irow,*icol,nrows,ncols; 1364 PetscScalar *av,*bv,*v = mat->v; 1365 Mat newmat; 1366 1367 PetscFunctionBegin; 1368 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1369 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1370 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1371 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1372 1373 /* Check submatrixcall */ 1374 if (scall == MAT_REUSE_MATRIX) { 1375 PetscInt n_cols,n_rows; 1376 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1377 if (n_rows != nrows || n_cols != ncols) { 1378 /* resize the result result matrix to match number of requested rows/columns */ 1379 ierr = MatSetSizes(*B,nrows,nrows,nrows,nrows);CHKERRQ(ierr); 1380 } 1381 newmat = *B; 1382 } else { 1383 /* Create and fill new matrix */ 1384 ierr = MatCreate(((PetscObject)A)->comm,&newmat);CHKERRQ(ierr); 1385 ierr = MatSetSizes(newmat,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1386 ierr = MatSetType(newmat,((PetscObject)A)->type_name);CHKERRQ(ierr); 1387 ierr = MatSeqDenseSetPreallocation(newmat,PETSC_NULL);CHKERRQ(ierr); 1388 } 1389 1390 /* Now extract the data pointers and do the copy,column at a time */ 1391 bv = ((Mat_SeqDense*)newmat->data)->v; 1392 1393 for (i=0; i<ncols; i++) { 1394 av = v + mat->lda*icol[i]; 1395 for (j=0; j<nrows; j++) { 1396 *bv++ = av[irow[j]]; 1397 } 1398 } 1399 1400 /* Assemble the matrices so that the correct flags are set */ 1401 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1402 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1403 1404 /* Free work space */ 1405 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1406 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1407 *B = newmat; 1408 PetscFunctionReturn(0); 1409 } 1410 1411 #undef __FUNCT__ 1412 #define __FUNCT__ "MatGetSubMatrices_SeqDense" 1413 PetscErrorCode MatGetSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1414 { 1415 PetscErrorCode ierr; 1416 PetscInt i; 1417 1418 PetscFunctionBegin; 1419 if (scall == MAT_INITIAL_MATRIX) { 1420 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1421 } 1422 1423 for (i=0; i<n; i++) { 1424 ierr = MatGetSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1425 } 1426 PetscFunctionReturn(0); 1427 } 1428 1429 #undef __FUNCT__ 1430 #define __FUNCT__ "MatAssemblyBegin_SeqDense" 1431 PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode) 1432 { 1433 PetscFunctionBegin; 1434 PetscFunctionReturn(0); 1435 } 1436 1437 #undef __FUNCT__ 1438 #define __FUNCT__ "MatAssemblyEnd_SeqDense" 1439 PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode) 1440 { 1441 PetscFunctionBegin; 1442 PetscFunctionReturn(0); 1443 } 1444 1445 #undef __FUNCT__ 1446 #define __FUNCT__ "MatCopy_SeqDense" 1447 PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str) 1448 { 1449 Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense *)B->data; 1450 PetscErrorCode ierr; 1451 PetscInt lda1=a->lda,lda2=b->lda, m=A->rmap.n,n=A->cmap.n, j; 1452 1453 PetscFunctionBegin; 1454 /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */ 1455 if (A->ops->copy != B->ops->copy) { 1456 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1457 PetscFunctionReturn(0); 1458 } 1459 if (m != B->rmap.n || n != B->cmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"size(B) != size(A)"); 1460 if (lda1>m || lda2>m) { 1461 for (j=0; j<n; j++) { 1462 ierr = PetscMemcpy(b->v+j*lda2,a->v+j*lda1,m*sizeof(PetscScalar));CHKERRQ(ierr); 1463 } 1464 } else { 1465 ierr = PetscMemcpy(b->v,a->v,A->rmap.n*A->cmap.n*sizeof(PetscScalar));CHKERRQ(ierr); 1466 } 1467 PetscFunctionReturn(0); 1468 } 1469 1470 #undef __FUNCT__ 1471 #define __FUNCT__ "MatSetUpPreallocation_SeqDense" 1472 PetscErrorCode MatSetUpPreallocation_SeqDense(Mat A) 1473 { 1474 PetscErrorCode ierr; 1475 1476 PetscFunctionBegin; 1477 ierr = MatSeqDenseSetPreallocation(A,0);CHKERRQ(ierr); 1478 PetscFunctionReturn(0); 1479 } 1480 1481 #undef __FUNCT__ 1482 #define __FUNCT__ "MatSetSizes_SeqDense" 1483 PetscErrorCode MatSetSizes_SeqDense(Mat A,PetscInt m,PetscInt n,PetscInt M,PetscInt N) 1484 { 1485 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1486 PetscErrorCode ierr; 1487 PetscFunctionBegin; 1488 /* this will not be called before lda, Mmax, and Nmax have been set */ 1489 m = PetscMax(m,M); 1490 n = PetscMax(n,N); 1491 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); 1492 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); 1493 A->rmap.n = A->rmap.n = m; 1494 A->cmap.n = A->cmap.N = n; 1495 if (a->changelda) a->lda = m; 1496 ierr = PetscMemzero(a->v,m*n*sizeof(PetscScalar));CHKERRQ(ierr); 1497 PetscFunctionReturn(0); 1498 } 1499 1500 /* ----------------------------------------------------------------*/ 1501 #undef __FUNCT__ 1502 #define __FUNCT__ "MatMatMult_SeqDense_SeqDense" 1503 PetscErrorCode MatMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 1504 { 1505 PetscErrorCode ierr; 1506 1507 PetscFunctionBegin; 1508 if (scall == MAT_INITIAL_MATRIX){ 1509 ierr = MatMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 1510 } 1511 ierr = MatMatMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 1512 PetscFunctionReturn(0); 1513 } 1514 1515 #undef __FUNCT__ 1516 #define __FUNCT__ "MatMatMultSymbolic_SeqDense_SeqDense" 1517 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 1518 { 1519 PetscErrorCode ierr; 1520 PetscInt m=A->rmap.n,n=B->cmap.n; 1521 Mat Cmat; 1522 1523 PetscFunctionBegin; 1524 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); 1525 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 1526 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 1527 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 1528 ierr = MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);CHKERRQ(ierr); 1529 Cmat->assembled = PETSC_TRUE; 1530 *C = Cmat; 1531 PetscFunctionReturn(0); 1532 } 1533 1534 #undef __FUNCT__ 1535 #define __FUNCT__ "MatMatMultNumeric_SeqDense_SeqDense" 1536 PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 1537 { 1538 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1539 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1540 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 1541 PetscBLASInt m,n,k; 1542 PetscScalar _DOne=1.0,_DZero=0.0; 1543 1544 PetscFunctionBegin; 1545 m = PetscBLASIntCast(A->rmap.n); 1546 n = PetscBLASIntCast(B->cmap.n); 1547 k = PetscBLASIntCast(A->cmap.n); 1548 BLASgemm_("N","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda); 1549 PetscFunctionReturn(0); 1550 } 1551 1552 #undef __FUNCT__ 1553 #define __FUNCT__ "MatMatMultTranspose_SeqDense_SeqDense" 1554 PetscErrorCode MatMatMultTranspose_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 1555 { 1556 PetscErrorCode ierr; 1557 1558 PetscFunctionBegin; 1559 if (scall == MAT_INITIAL_MATRIX){ 1560 ierr = MatMatMultTransposeSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 1561 } 1562 ierr = MatMatMultTransposeNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 1563 PetscFunctionReturn(0); 1564 } 1565 1566 #undef __FUNCT__ 1567 #define __FUNCT__ "MatMatMultTransposeSymbolic_SeqDense_SeqDense" 1568 PetscErrorCode MatMatMultTransposeSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 1569 { 1570 PetscErrorCode ierr; 1571 PetscInt m=A->cmap.n,n=B->cmap.n; 1572 Mat Cmat; 1573 1574 PetscFunctionBegin; 1575 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); 1576 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 1577 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 1578 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 1579 ierr = MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);CHKERRQ(ierr); 1580 Cmat->assembled = PETSC_TRUE; 1581 *C = Cmat; 1582 PetscFunctionReturn(0); 1583 } 1584 1585 #undef __FUNCT__ 1586 #define __FUNCT__ "MatMatMultTransposeNumeric_SeqDense_SeqDense" 1587 PetscErrorCode MatMatMultTransposeNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 1588 { 1589 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1590 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1591 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 1592 PetscBLASInt m,n,k; 1593 PetscScalar _DOne=1.0,_DZero=0.0; 1594 1595 PetscFunctionBegin; 1596 m = PetscBLASIntCast(A->rmap.n); 1597 n = PetscBLASIntCast(B->cmap.n); 1598 k = PetscBLASIntCast(A->cmap.n); 1599 BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda); 1600 PetscFunctionReturn(0); 1601 } 1602 1603 #undef __FUNCT__ 1604 #define __FUNCT__ "MatGetRowMax_SeqDense" 1605 PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[]) 1606 { 1607 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1608 PetscErrorCode ierr; 1609 PetscInt i,j,m = A->rmap.n,n = A->cmap.n,p; 1610 PetscScalar *x; 1611 MatScalar *aa = a->v; 1612 1613 PetscFunctionBegin; 1614 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1615 1616 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1617 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1618 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1619 if (p != A->rmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1620 for (i=0; i<m; i++) { 1621 x[i] = aa[i]; if (idx) idx[i] = 0; 1622 for (j=1; j<n; j++){ 1623 if (PetscRealPart(x[i]) < PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;} 1624 } 1625 } 1626 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1627 PetscFunctionReturn(0); 1628 } 1629 1630 #undef __FUNCT__ 1631 #define __FUNCT__ "MatGetRowMaxAbs_SeqDense" 1632 PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[]) 1633 { 1634 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1635 PetscErrorCode ierr; 1636 PetscInt i,j,m = A->rmap.n,n = A->cmap.n,p; 1637 PetscScalar *x; 1638 PetscReal atmp; 1639 MatScalar *aa = a->v; 1640 1641 PetscFunctionBegin; 1642 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1643 1644 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1645 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1646 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1647 if (p != A->rmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1648 for (i=0; i<m; i++) { 1649 x[i] = PetscAbsScalar(aa[i]); if (idx) idx[i] = 0; 1650 for (j=1; j<n; j++){ 1651 atmp = PetscAbsScalar(aa[i+m*j]); 1652 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;} 1653 } 1654 } 1655 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1656 PetscFunctionReturn(0); 1657 } 1658 1659 #undef __FUNCT__ 1660 #define __FUNCT__ "MatGetRowMin_SeqDense" 1661 PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[]) 1662 { 1663 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1664 PetscErrorCode ierr; 1665 PetscInt i,j,m = A->rmap.n,n = A->cmap.n,p; 1666 PetscScalar *x; 1667 MatScalar *aa = a->v; 1668 1669 PetscFunctionBegin; 1670 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1671 1672 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1673 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1674 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1675 if (p != A->rmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1676 for (i=0; i<m; i++) { 1677 x[i] = aa[i]; if (idx) idx[i] = 0; 1678 for (j=1; j<n; j++){ 1679 if (PetscRealPart(x[i]) > PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;} 1680 } 1681 } 1682 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1683 PetscFunctionReturn(0); 1684 } 1685 1686 #undef __FUNCT__ 1687 #define __FUNCT__ "MatGetColumnVector_SeqDense" 1688 PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col) 1689 { 1690 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1691 PetscErrorCode ierr; 1692 PetscScalar *x; 1693 1694 PetscFunctionBegin; 1695 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1696 1697 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1698 ierr = PetscMemcpy(x,a->v+col*a->lda,A->rmap.n*sizeof(PetscScalar));CHKERRQ(ierr); 1699 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1700 PetscFunctionReturn(0); 1701 } 1702 1703 /* -------------------------------------------------------------------*/ 1704 static struct _MatOps MatOps_Values = {MatSetValues_SeqDense, 1705 MatGetRow_SeqDense, 1706 MatRestoreRow_SeqDense, 1707 MatMult_SeqDense, 1708 /* 4*/ MatMultAdd_SeqDense, 1709 MatMultTranspose_SeqDense, 1710 MatMultTransposeAdd_SeqDense, 1711 MatSolve_SeqDense, 1712 MatSolveAdd_SeqDense, 1713 MatSolveTranspose_SeqDense, 1714 /*10*/ MatSolveTransposeAdd_SeqDense, 1715 MatLUFactor_SeqDense, 1716 MatCholeskyFactor_SeqDense, 1717 MatRelax_SeqDense, 1718 MatTranspose_SeqDense, 1719 /*15*/ MatGetInfo_SeqDense, 1720 MatEqual_SeqDense, 1721 MatGetDiagonal_SeqDense, 1722 MatDiagonalScale_SeqDense, 1723 MatNorm_SeqDense, 1724 /*20*/ MatAssemblyBegin_SeqDense, 1725 MatAssemblyEnd_SeqDense, 1726 0, 1727 MatSetOption_SeqDense, 1728 MatZeroEntries_SeqDense, 1729 /*25*/ MatZeroRows_SeqDense, 1730 MatLUFactorSymbolic_SeqDense, 1731 MatLUFactorNumeric_SeqDense, 1732 MatCholeskyFactorSymbolic_SeqDense, 1733 MatCholeskyFactorNumeric_SeqDense, 1734 /*30*/ MatSetUpPreallocation_SeqDense, 1735 0, 1736 0, 1737 MatGetArray_SeqDense, 1738 MatRestoreArray_SeqDense, 1739 /*35*/ MatDuplicate_SeqDense, 1740 0, 1741 0, 1742 0, 1743 0, 1744 /*40*/ MatAXPY_SeqDense, 1745 MatGetSubMatrices_SeqDense, 1746 0, 1747 MatGetValues_SeqDense, 1748 MatCopy_SeqDense, 1749 /*45*/ MatGetRowMax_SeqDense, 1750 MatScale_SeqDense, 1751 0, 1752 0, 1753 0, 1754 /*50*/ 0, 1755 0, 1756 0, 1757 0, 1758 0, 1759 /*55*/ 0, 1760 0, 1761 0, 1762 0, 1763 0, 1764 /*60*/ 0, 1765 MatDestroy_SeqDense, 1766 MatView_SeqDense, 1767 0, 1768 0, 1769 /*65*/ 0, 1770 0, 1771 0, 1772 0, 1773 0, 1774 /*70*/ MatGetRowMaxAbs_SeqDense, 1775 0, 1776 0, 1777 0, 1778 0, 1779 /*75*/ 0, 1780 0, 1781 0, 1782 0, 1783 0, 1784 /*80*/ 0, 1785 0, 1786 0, 1787 0, 1788 /*84*/ MatLoad_SeqDense, 1789 0, 1790 MatIsHermitian_SeqDense, 1791 0, 1792 0, 1793 0, 1794 /*90*/ MatMatMult_SeqDense_SeqDense, 1795 MatMatMultSymbolic_SeqDense_SeqDense, 1796 MatMatMultNumeric_SeqDense_SeqDense, 1797 0, 1798 0, 1799 /*95*/ 0, 1800 MatMatMultTranspose_SeqDense_SeqDense, 1801 MatMatMultTransposeSymbolic_SeqDense_SeqDense, 1802 MatMatMultTransposeNumeric_SeqDense_SeqDense, 1803 0, 1804 /*100*/0, 1805 0, 1806 0, 1807 0, 1808 MatSetSizes_SeqDense, 1809 0, 1810 0, 1811 0, 1812 0, 1813 0, 1814 /*110*/0, 1815 0, 1816 MatGetRowMin_SeqDense, 1817 MatGetColumnVector_SeqDense 1818 }; 1819 1820 #undef __FUNCT__ 1821 #define __FUNCT__ "MatCreateSeqDense" 1822 /*@C 1823 MatCreateSeqDense - Creates a sequential dense matrix that 1824 is stored in column major order (the usual Fortran 77 manner). Many 1825 of the matrix operations use the BLAS and LAPACK routines. 1826 1827 Collective on MPI_Comm 1828 1829 Input Parameters: 1830 + comm - MPI communicator, set to PETSC_COMM_SELF 1831 . m - number of rows 1832 . n - number of columns 1833 - data - optional location of matrix data. Set data=PETSC_NULL for PETSc 1834 to control all matrix memory allocation. 1835 1836 Output Parameter: 1837 . A - the matrix 1838 1839 Notes: 1840 The data input variable is intended primarily for Fortran programmers 1841 who wish to allocate their own matrix memory space. Most users should 1842 set data=PETSC_NULL. 1843 1844 Level: intermediate 1845 1846 .keywords: dense, matrix, LAPACK, BLAS 1847 1848 .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues() 1849 @*/ 1850 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A) 1851 { 1852 PetscErrorCode ierr; 1853 1854 PetscFunctionBegin; 1855 ierr = MatCreate(comm,A);CHKERRQ(ierr); 1856 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 1857 ierr = MatSetType(*A,MATSEQDENSE);CHKERRQ(ierr); 1858 ierr = MatSeqDenseSetPreallocation(*A,data);CHKERRQ(ierr); 1859 PetscFunctionReturn(0); 1860 } 1861 1862 #undef __FUNCT__ 1863 #define __FUNCT__ "MatSeqDenseSetPreallocation" 1864 /*@C 1865 MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements 1866 1867 Collective on MPI_Comm 1868 1869 Input Parameters: 1870 + A - the matrix 1871 - data - the array (or PETSC_NULL) 1872 1873 Notes: 1874 The data input variable is intended primarily for Fortran programmers 1875 who wish to allocate their own matrix memory space. Most users should 1876 need not call this routine. 1877 1878 Level: intermediate 1879 1880 .keywords: dense, matrix, LAPACK, BLAS 1881 1882 .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues() 1883 @*/ 1884 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqDenseSetPreallocation(Mat B,PetscScalar data[]) 1885 { 1886 PetscErrorCode ierr,(*f)(Mat,PetscScalar[]); 1887 1888 PetscFunctionBegin; 1889 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 1890 if (f) { 1891 ierr = (*f)(B,data);CHKERRQ(ierr); 1892 } 1893 PetscFunctionReturn(0); 1894 } 1895 1896 EXTERN_C_BEGIN 1897 #undef __FUNCT__ 1898 #define __FUNCT__ "MatSeqDenseSetPreallocation_SeqDense" 1899 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data) 1900 { 1901 Mat_SeqDense *b; 1902 PetscErrorCode ierr; 1903 1904 PetscFunctionBegin; 1905 B->preallocated = PETSC_TRUE; 1906 b = (Mat_SeqDense*)B->data; 1907 if (!data) { /* petsc-allocated storage */ 1908 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 1909 ierr = PetscMalloc((b->lda*b->Nmax+1)*sizeof(PetscScalar),&b->v);CHKERRQ(ierr); 1910 ierr = PetscMemzero(b->v,b->lda*b->Nmax*sizeof(PetscScalar));CHKERRQ(ierr); 1911 ierr = PetscLogObjectMemory(B,b->lda*b->Nmax*sizeof(PetscScalar));CHKERRQ(ierr); 1912 b->user_alloc = PETSC_FALSE; 1913 } else { /* user-allocated storage */ 1914 if (!b->user_alloc) { ierr = PetscFree(b->v);CHKERRQ(ierr); } 1915 b->v = data; 1916 b->user_alloc = PETSC_TRUE; 1917 } 1918 PetscFunctionReturn(0); 1919 } 1920 EXTERN_C_END 1921 1922 #undef __FUNCT__ 1923 #define __FUNCT__ "MatSeqDenseSetLDA" 1924 /*@C 1925 MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array 1926 1927 Input parameter: 1928 + A - the matrix 1929 - lda - the leading dimension 1930 1931 Notes: 1932 This routine is to be used in conjunction with MatSeqDenseSetPreallocation; 1933 it asserts that the preallocation has a leading dimension (the LDA parameter 1934 of Blas and Lapack fame) larger than M, the first dimension of the matrix. 1935 1936 Level: intermediate 1937 1938 .keywords: dense, matrix, LAPACK, BLAS 1939 1940 .seealso: MatCreate(), MatCreateSeqDense(), MatSeqDenseSetPreallocation(), MatSetMaximumSize() 1941 @*/ 1942 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqDenseSetLDA(Mat B,PetscInt lda) 1943 { 1944 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1945 1946 PetscFunctionBegin; 1947 if (lda < B->rmap.n) SETERRQ2(PETSC_ERR_ARG_SIZ,"LDA %D must be at least matrix dimension %D",lda,B->rmap.n); 1948 b->lda = lda; 1949 b->changelda = PETSC_FALSE; 1950 b->Mmax = PetscMax(b->Mmax,lda); 1951 PetscFunctionReturn(0); 1952 } 1953 1954 /*MC 1955 MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices. 1956 1957 Options Database Keys: 1958 . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions() 1959 1960 Level: beginner 1961 1962 .seealso: MatCreateSeqDense() 1963 1964 M*/ 1965 1966 EXTERN_C_BEGIN 1967 #undef __FUNCT__ 1968 #define __FUNCT__ "MatCreate_SeqDense" 1969 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_SeqDense(Mat B) 1970 { 1971 Mat_SeqDense *b; 1972 PetscErrorCode ierr; 1973 PetscMPIInt size; 1974 1975 PetscFunctionBegin; 1976 ierr = MPI_Comm_size(((PetscObject)B)->comm,&size);CHKERRQ(ierr); 1977 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 1978 1979 B->rmap.bs = B->cmap.bs = 1; 1980 ierr = PetscMapSetUp(&B->rmap);CHKERRQ(ierr); 1981 ierr = PetscMapSetUp(&B->cmap);CHKERRQ(ierr); 1982 1983 ierr = PetscNewLog(B,Mat_SeqDense,&b);CHKERRQ(ierr); 1984 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1985 B->factor = 0; 1986 B->mapping = 0; 1987 B->data = (void*)b; 1988 1989 1990 b->pivots = 0; 1991 b->roworiented = PETSC_TRUE; 1992 b->v = 0; 1993 b->lda = B->rmap.n; 1994 b->changelda = PETSC_FALSE; 1995 b->Mmax = B->rmap.n; 1996 b->Nmax = B->cmap.n; 1997 1998 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqDenseSetPreallocation_C", 1999 "MatSeqDenseSetPreallocation_SeqDense", 2000 MatSeqDenseSetPreallocation_SeqDense);CHKERRQ(ierr); 2001 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMult_seqaij_seqdense_C", 2002 "MatMatMult_SeqAIJ_SeqDense", 2003 MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 2004 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C", 2005 "MatMatMultSymbolic_SeqAIJ_SeqDense", 2006 MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 2007 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C", 2008 "MatMatMultNumeric_SeqAIJ_SeqDense", 2009 MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 2010 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);CHKERRQ(ierr); 2011 PetscFunctionReturn(0); 2012 } 2013 2014 2015 EXTERN_C_END 2016