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