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 667 PetscFunctionBegin; 668 /* row-oriented output */ 669 for (i=0; i<m; i++) { 670 if (indexm[i] < 0) {v += n;continue;} 671 for (j=0; j<n; j++) { 672 if (indexn[i] < 0) {v++; continue;} 673 *v++ = mat->v[indexn[j]*mat->lda + indexm[i]]; 674 } 675 } 676 PetscFunctionReturn(0); 677 } 678 679 /* -----------------------------------------------------------------*/ 680 681 #include "petscsys.h" 682 683 #undef __FUNCT__ 684 #define __FUNCT__ "MatLoad_SeqDense" 685 PetscErrorCode MatLoad_SeqDense(PetscViewer viewer, MatType type,Mat *A) 686 { 687 Mat_SeqDense *a; 688 Mat B; 689 PetscErrorCode ierr; 690 PetscInt *scols,i,j,nz,header[4]; 691 int fd; 692 PetscMPIInt size; 693 PetscInt *rowlengths = 0,M,N,*cols; 694 PetscScalar *vals,*svals,*v,*w; 695 MPI_Comm comm = ((PetscObject)viewer)->comm; 696 697 PetscFunctionBegin; 698 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 699 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"view must have one processor"); 700 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 701 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 702 if (header[0] != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Not matrix object"); 703 M = header[1]; N = header[2]; nz = header[3]; 704 705 if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense */ 706 ierr = MatCreate(comm,A);CHKERRQ(ierr); 707 ierr = MatSetSizes(*A,M,N,M,N);CHKERRQ(ierr); 708 ierr = MatSetType(*A,type);CHKERRQ(ierr); 709 ierr = MatSeqDenseSetPreallocation(*A,PETSC_NULL);CHKERRQ(ierr); 710 B = *A; 711 a = (Mat_SeqDense*)B->data; 712 v = a->v; 713 /* Allocate some temp space to read in the values and then flip them 714 from row major to column major */ 715 ierr = PetscMalloc((M*N > 0 ? M*N : 1)*sizeof(PetscScalar),&w);CHKERRQ(ierr); 716 /* read in nonzero values */ 717 ierr = PetscBinaryRead(fd,w,M*N,PETSC_SCALAR);CHKERRQ(ierr); 718 /* now flip the values and store them in the matrix*/ 719 for (j=0; j<N; j++) { 720 for (i=0; i<M; i++) { 721 *v++ =w[i*N+j]; 722 } 723 } 724 ierr = PetscFree(w);CHKERRQ(ierr); 725 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 726 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 727 } else { 728 /* read row lengths */ 729 ierr = PetscMalloc((M+1)*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 730 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 731 732 /* create our matrix */ 733 ierr = MatCreate(comm,A);CHKERRQ(ierr); 734 ierr = MatSetSizes(*A,M,N,M,N);CHKERRQ(ierr); 735 ierr = MatSetType(*A,type);CHKERRQ(ierr); 736 ierr = MatSeqDenseSetPreallocation(*A,PETSC_NULL);CHKERRQ(ierr); 737 B = *A; 738 a = (Mat_SeqDense*)B->data; 739 v = a->v; 740 741 /* read column indices and nonzeros */ 742 ierr = PetscMalloc((nz+1)*sizeof(PetscInt),&scols);CHKERRQ(ierr); 743 cols = scols; 744 ierr = PetscBinaryRead(fd,cols,nz,PETSC_INT);CHKERRQ(ierr); 745 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&svals);CHKERRQ(ierr); 746 vals = svals; 747 ierr = PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);CHKERRQ(ierr); 748 749 /* insert into matrix */ 750 for (i=0; i<M; i++) { 751 for (j=0; j<rowlengths[i]; j++) v[i+M*scols[j]] = svals[j]; 752 svals += rowlengths[i]; scols += rowlengths[i]; 753 } 754 ierr = PetscFree(vals);CHKERRQ(ierr); 755 ierr = PetscFree(cols);CHKERRQ(ierr); 756 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 757 758 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 759 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 760 } 761 PetscFunctionReturn(0); 762 } 763 764 #include "petscsys.h" 765 766 #undef __FUNCT__ 767 #define __FUNCT__ "MatView_SeqDense_ASCII" 768 static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer) 769 { 770 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 771 PetscErrorCode ierr; 772 PetscInt i,j; 773 const char *name; 774 PetscScalar *v; 775 PetscViewerFormat format; 776 #if defined(PETSC_USE_COMPLEX) 777 PetscTruth allreal = PETSC_TRUE; 778 #endif 779 780 PetscFunctionBegin; 781 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 782 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 783 if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 784 PetscFunctionReturn(0); /* do nothing for now */ 785 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 786 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 787 for (i=0; i<A->rmap.n; i++) { 788 v = a->v + i; 789 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 790 for (j=0; j<A->cmap.n; j++) { 791 #if defined(PETSC_USE_COMPLEX) 792 if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) { 793 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i) ",j,PetscRealPart(*v),PetscImaginaryPart(*v));CHKERRQ(ierr); 794 } else if (PetscRealPart(*v)) { 795 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",j,PetscRealPart(*v));CHKERRQ(ierr); 796 } 797 #else 798 if (*v) { 799 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",j,*v);CHKERRQ(ierr); 800 } 801 #endif 802 v += a->lda; 803 } 804 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 805 } 806 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 807 } else { 808 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 809 #if defined(PETSC_USE_COMPLEX) 810 /* determine if matrix has all real values */ 811 v = a->v; 812 for (i=0; i<A->rmap.n*A->cmap.n; i++) { 813 if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break ;} 814 } 815 #endif 816 if (format == PETSC_VIEWER_ASCII_MATLAB) { 817 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 818 ierr = PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",A->rmap.n,A->cmap.n);CHKERRQ(ierr); 819 ierr = PetscViewerASCIIPrintf(viewer,"%s = zeros(%D,%D);\n",name,A->rmap.n,A->cmap.n);CHKERRQ(ierr); 820 ierr = PetscViewerASCIIPrintf(viewer,"%s = [\n",name);CHKERRQ(ierr); 821 } 822 823 for (i=0; i<A->rmap.n; i++) { 824 v = a->v + i; 825 for (j=0; j<A->cmap.n; j++) { 826 #if defined(PETSC_USE_COMPLEX) 827 if (allreal) { 828 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e ",PetscRealPart(*v));CHKERRQ(ierr); 829 } else { 830 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e + %18.16e i ",PetscRealPart(*v),PetscImaginaryPart(*v));CHKERRQ(ierr); 831 } 832 #else 833 ierr = PetscViewerASCIIPrintf(viewer,"%18.16e ",*v);CHKERRQ(ierr); 834 #endif 835 v += a->lda; 836 } 837 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 838 } 839 if (format == PETSC_VIEWER_ASCII_MATLAB) { 840 ierr = PetscViewerASCIIPrintf(viewer,"];\n");CHKERRQ(ierr); 841 } 842 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 843 } 844 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 845 PetscFunctionReturn(0); 846 } 847 848 #undef __FUNCT__ 849 #define __FUNCT__ "MatView_SeqDense_Binary" 850 static PetscErrorCode MatView_SeqDense_Binary(Mat A,PetscViewer viewer) 851 { 852 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 853 PetscErrorCode ierr; 854 int fd; 855 PetscInt ict,j,n = A->cmap.n,m = A->rmap.n,i,*col_lens,nz = m*n; 856 PetscScalar *v,*anonz,*vals; 857 PetscViewerFormat format; 858 859 PetscFunctionBegin; 860 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 861 862 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 863 if (format == PETSC_VIEWER_BINARY_NATIVE) { 864 /* store the matrix as a dense matrix */ 865 ierr = PetscMalloc(4*sizeof(PetscInt),&col_lens);CHKERRQ(ierr); 866 col_lens[0] = MAT_FILE_COOKIE; 867 col_lens[1] = m; 868 col_lens[2] = n; 869 col_lens[3] = MATRIX_BINARY_FORMAT_DENSE; 870 ierr = PetscBinaryWrite(fd,col_lens,4,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 871 ierr = PetscFree(col_lens);CHKERRQ(ierr); 872 873 /* write out matrix, by rows */ 874 ierr = PetscMalloc((m*n+1)*sizeof(PetscScalar),&vals);CHKERRQ(ierr); 875 v = a->v; 876 for (i=0; i<m; i++) { 877 for (j=0; j<n; j++) { 878 vals[i + j*m] = *v++; 879 } 880 } 881 ierr = PetscBinaryWrite(fd,vals,n*m,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 882 ierr = PetscFree(vals);CHKERRQ(ierr); 883 } else { 884 ierr = PetscMalloc((4+nz)*sizeof(PetscInt),&col_lens);CHKERRQ(ierr); 885 col_lens[0] = MAT_FILE_COOKIE; 886 col_lens[1] = m; 887 col_lens[2] = n; 888 col_lens[3] = nz; 889 890 /* store lengths of each row and write (including header) to file */ 891 for (i=0; i<m; i++) col_lens[4+i] = n; 892 ierr = PetscBinaryWrite(fd,col_lens,4+m,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 893 894 /* Possibly should write in smaller increments, not whole matrix at once? */ 895 /* store column indices (zero start index) */ 896 ict = 0; 897 for (i=0; i<m; i++) { 898 for (j=0; j<n; j++) col_lens[ict++] = j; 899 } 900 ierr = PetscBinaryWrite(fd,col_lens,nz,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 901 ierr = PetscFree(col_lens);CHKERRQ(ierr); 902 903 /* store nonzero values */ 904 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&anonz);CHKERRQ(ierr); 905 ict = 0; 906 for (i=0; i<m; i++) { 907 v = a->v + i; 908 for (j=0; j<n; j++) { 909 anonz[ict++] = *v; v += a->lda; 910 } 911 } 912 ierr = PetscBinaryWrite(fd,anonz,nz,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 913 ierr = PetscFree(anonz);CHKERRQ(ierr); 914 } 915 PetscFunctionReturn(0); 916 } 917 918 #undef __FUNCT__ 919 #define __FUNCT__ "MatView_SeqDense_Draw_Zoom" 920 PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa) 921 { 922 Mat A = (Mat) Aa; 923 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 924 PetscErrorCode ierr; 925 PetscInt m = A->rmap.n,n = A->cmap.n,color,i,j; 926 PetscScalar *v = a->v; 927 PetscViewer viewer; 928 PetscDraw popup; 929 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r,scale,maxv = 0.0; 930 PetscViewerFormat format; 931 932 PetscFunctionBegin; 933 934 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 935 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 936 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 937 938 /* Loop over matrix elements drawing boxes */ 939 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 940 /* Blue for negative and Red for positive */ 941 color = PETSC_DRAW_BLUE; 942 for(j = 0; j < n; j++) { 943 x_l = j; 944 x_r = x_l + 1.0; 945 for(i = 0; i < m; i++) { 946 y_l = m - i - 1.0; 947 y_r = y_l + 1.0; 948 #if defined(PETSC_USE_COMPLEX) 949 if (PetscRealPart(v[j*m+i]) > 0.) { 950 color = PETSC_DRAW_RED; 951 } else if (PetscRealPart(v[j*m+i]) < 0.) { 952 color = PETSC_DRAW_BLUE; 953 } else { 954 continue; 955 } 956 #else 957 if (v[j*m+i] > 0.) { 958 color = PETSC_DRAW_RED; 959 } else if (v[j*m+i] < 0.) { 960 color = PETSC_DRAW_BLUE; 961 } else { 962 continue; 963 } 964 #endif 965 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 966 } 967 } 968 } else { 969 /* use contour shading to indicate magnitude of values */ 970 /* first determine max of all nonzero values */ 971 for(i = 0; i < m*n; i++) { 972 if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]); 973 } 974 scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/maxv; 975 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 976 if (popup) {ierr = PetscDrawScalePopup(popup,0.0,maxv);CHKERRQ(ierr);} 977 for(j = 0; j < n; j++) { 978 x_l = j; 979 x_r = x_l + 1.0; 980 for(i = 0; i < m; i++) { 981 y_l = m - i - 1.0; 982 y_r = y_l + 1.0; 983 color = PETSC_DRAW_BASIC_COLORS + (int)(scale*PetscAbsScalar(v[j*m+i])); 984 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 985 } 986 } 987 } 988 PetscFunctionReturn(0); 989 } 990 991 #undef __FUNCT__ 992 #define __FUNCT__ "MatView_SeqDense_Draw" 993 PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer) 994 { 995 PetscDraw draw; 996 PetscTruth isnull; 997 PetscReal xr,yr,xl,yl,h,w; 998 PetscErrorCode ierr; 999 1000 PetscFunctionBegin; 1001 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 1002 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 1003 if (isnull) PetscFunctionReturn(0); 1004 1005 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 1006 xr = A->cmap.n; yr = A->rmap.n; h = yr/10.0; w = xr/10.0; 1007 xr += w; yr += h; xl = -w; yl = -h; 1008 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 1009 ierr = PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);CHKERRQ(ierr); 1010 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);CHKERRQ(ierr); 1011 PetscFunctionReturn(0); 1012 } 1013 1014 #undef __FUNCT__ 1015 #define __FUNCT__ "MatView_SeqDense" 1016 PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer) 1017 { 1018 PetscErrorCode ierr; 1019 PetscTruth iascii,isbinary,isdraw; 1020 1021 PetscFunctionBegin; 1022 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 1023 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_BINARY,&isbinary);CHKERRQ(ierr); 1024 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_DRAW,&isdraw);CHKERRQ(ierr); 1025 1026 if (iascii) { 1027 ierr = MatView_SeqDense_ASCII(A,viewer);CHKERRQ(ierr); 1028 } else if (isbinary) { 1029 ierr = MatView_SeqDense_Binary(A,viewer);CHKERRQ(ierr); 1030 } else if (isdraw) { 1031 ierr = MatView_SeqDense_Draw(A,viewer);CHKERRQ(ierr); 1032 } else { 1033 SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported by dense matrix",((PetscObject)viewer)->type_name); 1034 } 1035 PetscFunctionReturn(0); 1036 } 1037 1038 #undef __FUNCT__ 1039 #define __FUNCT__ "MatDestroy_SeqDense" 1040 PetscErrorCode MatDestroy_SeqDense(Mat mat) 1041 { 1042 Mat_SeqDense *l = (Mat_SeqDense*)mat->data; 1043 PetscErrorCode ierr; 1044 1045 PetscFunctionBegin; 1046 #if defined(PETSC_USE_LOG) 1047 PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->rmap.n,mat->cmap.n); 1048 #endif 1049 ierr = PetscFree(l->pivots);CHKERRQ(ierr); 1050 if (!l->user_alloc) {ierr = PetscFree(l->v);CHKERRQ(ierr);} 1051 ierr = PetscFree(l);CHKERRQ(ierr); 1052 1053 ierr = PetscObjectChangeTypeName((PetscObject)mat,0);CHKERRQ(ierr); 1054 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatSeqDenseSetPreallocation_C","",PETSC_NULL);CHKERRQ(ierr); 1055 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatMatMult_seqaij_seqdense_C","",PETSC_NULL);CHKERRQ(ierr); 1056 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatMatMultSymbolic_seqaij_seqdense_C","",PETSC_NULL);CHKERRQ(ierr); 1057 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatMatMultNumeric_seqaij_seqdense_C","",PETSC_NULL);CHKERRQ(ierr); 1058 PetscFunctionReturn(0); 1059 } 1060 1061 #undef __FUNCT__ 1062 #define __FUNCT__ "MatTranspose_SeqDense" 1063 PetscErrorCode MatTranspose_SeqDense(Mat A,Mat *matout) 1064 { 1065 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1066 PetscErrorCode ierr; 1067 PetscInt k,j,m,n,M; 1068 PetscScalar *v,tmp; 1069 1070 PetscFunctionBegin; 1071 v = mat->v; m = A->rmap.n; M = mat->lda; n = A->cmap.n; 1072 if (!matout) { /* in place transpose */ 1073 if (m != n) { 1074 SETERRQ(PETSC_ERR_SUP,"Can not transpose non-square matrix in place"); 1075 } else { 1076 for (j=0; j<m; j++) { 1077 for (k=0; k<j; k++) { 1078 tmp = v[j + k*M]; 1079 v[j + k*M] = v[k + j*M]; 1080 v[k + j*M] = tmp; 1081 } 1082 } 1083 } 1084 } else { /* out-of-place transpose */ 1085 Mat tmat; 1086 Mat_SeqDense *tmatd; 1087 PetscScalar *v2; 1088 1089 ierr = MatCreate(A->comm,&tmat);CHKERRQ(ierr); 1090 ierr = MatSetSizes(tmat,A->cmap.n,A->rmap.n,A->cmap.n,A->rmap.n);CHKERRQ(ierr); 1091 ierr = MatSetType(tmat,A->type_name);CHKERRQ(ierr); 1092 ierr = MatSeqDenseSetPreallocation(tmat,PETSC_NULL);CHKERRQ(ierr); 1093 tmatd = (Mat_SeqDense*)tmat->data; 1094 v = mat->v; v2 = tmatd->v; 1095 for (j=0; j<n; j++) { 1096 for (k=0; k<m; k++) v2[j + k*n] = v[k + j*M]; 1097 } 1098 ierr = MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1099 ierr = MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1100 *matout = tmat; 1101 } 1102 PetscFunctionReturn(0); 1103 } 1104 1105 #undef __FUNCT__ 1106 #define __FUNCT__ "MatEqual_SeqDense" 1107 PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscTruth *flg) 1108 { 1109 Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data; 1110 Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data; 1111 PetscInt i,j; 1112 PetscScalar *v1 = mat1->v,*v2 = mat2->v; 1113 1114 PetscFunctionBegin; 1115 if (A1->rmap.n != A2->rmap.n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1116 if (A1->cmap.n != A2->cmap.n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1117 for (i=0; i<A1->rmap.n; i++) { 1118 v1 = mat1->v+i; v2 = mat2->v+i; 1119 for (j=0; j<A1->cmap.n; j++) { 1120 if (*v1 != *v2) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1121 v1 += mat1->lda; v2 += mat2->lda; 1122 } 1123 } 1124 *flg = PETSC_TRUE; 1125 PetscFunctionReturn(0); 1126 } 1127 1128 #undef __FUNCT__ 1129 #define __FUNCT__ "MatGetDiagonal_SeqDense" 1130 PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v) 1131 { 1132 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1133 PetscErrorCode ierr; 1134 PetscInt i,n,len; 1135 PetscScalar *x,zero = 0.0; 1136 1137 PetscFunctionBegin; 1138 ierr = VecSet(v,zero);CHKERRQ(ierr); 1139 ierr = VecGetSize(v,&n);CHKERRQ(ierr); 1140 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1141 len = PetscMin(A->rmap.n,A->cmap.n); 1142 if (n != A->rmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec"); 1143 for (i=0; i<len; i++) { 1144 x[i] = mat->v[i*mat->lda + i]; 1145 } 1146 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1147 PetscFunctionReturn(0); 1148 } 1149 1150 #undef __FUNCT__ 1151 #define __FUNCT__ "MatDiagonalScale_SeqDense" 1152 PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr) 1153 { 1154 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1155 PetscScalar *l,*r,x,*v; 1156 PetscErrorCode ierr; 1157 PetscInt i,j,m = A->rmap.n,n = A->cmap.n; 1158 1159 PetscFunctionBegin; 1160 if (ll) { 1161 ierr = VecGetSize(ll,&m);CHKERRQ(ierr); 1162 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1163 if (m != A->rmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size"); 1164 for (i=0; i<m; i++) { 1165 x = l[i]; 1166 v = mat->v + i; 1167 for (j=0; j<n; j++) { (*v) *= x; v+= m;} 1168 } 1169 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1170 ierr = PetscLogFlops(n*m);CHKERRQ(ierr); 1171 } 1172 if (rr) { 1173 ierr = VecGetSize(rr,&n);CHKERRQ(ierr); 1174 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1175 if (n != A->cmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size"); 1176 for (i=0; i<n; i++) { 1177 x = r[i]; 1178 v = mat->v + i*m; 1179 for (j=0; j<m; j++) { (*v++) *= x;} 1180 } 1181 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1182 ierr = PetscLogFlops(n*m);CHKERRQ(ierr); 1183 } 1184 PetscFunctionReturn(0); 1185 } 1186 1187 #undef __FUNCT__ 1188 #define __FUNCT__ "MatNorm_SeqDense" 1189 PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm) 1190 { 1191 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1192 PetscScalar *v = mat->v; 1193 PetscReal sum = 0.0; 1194 PetscInt lda=mat->lda,m=A->rmap.n,i,j; 1195 PetscErrorCode ierr; 1196 1197 PetscFunctionBegin; 1198 if (type == NORM_FROBENIUS) { 1199 if (lda>m) { 1200 for (j=0; j<A->cmap.n; j++) { 1201 v = mat->v+j*lda; 1202 for (i=0; i<m; i++) { 1203 #if defined(PETSC_USE_COMPLEX) 1204 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1205 #else 1206 sum += (*v)*(*v); v++; 1207 #endif 1208 } 1209 } 1210 } else { 1211 for (i=0; i<A->cmap.n*A->rmap.n; i++) { 1212 #if defined(PETSC_USE_COMPLEX) 1213 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1214 #else 1215 sum += (*v)*(*v); v++; 1216 #endif 1217 } 1218 } 1219 *nrm = sqrt(sum); 1220 ierr = PetscLogFlops(2*A->cmap.n*A->rmap.n);CHKERRQ(ierr); 1221 } else if (type == NORM_1) { 1222 *nrm = 0.0; 1223 for (j=0; j<A->cmap.n; j++) { 1224 v = mat->v + j*mat->lda; 1225 sum = 0.0; 1226 for (i=0; i<A->rmap.n; i++) { 1227 sum += PetscAbsScalar(*v); v++; 1228 } 1229 if (sum > *nrm) *nrm = sum; 1230 } 1231 ierr = PetscLogFlops(A->cmap.n*A->rmap.n);CHKERRQ(ierr); 1232 } else if (type == NORM_INFINITY) { 1233 *nrm = 0.0; 1234 for (j=0; j<A->rmap.n; j++) { 1235 v = mat->v + j; 1236 sum = 0.0; 1237 for (i=0; i<A->cmap.n; i++) { 1238 sum += PetscAbsScalar(*v); v += mat->lda; 1239 } 1240 if (sum > *nrm) *nrm = sum; 1241 } 1242 ierr = PetscLogFlops(A->cmap.n*A->rmap.n);CHKERRQ(ierr); 1243 } else { 1244 SETERRQ(PETSC_ERR_SUP,"No two norm"); 1245 } 1246 PetscFunctionReturn(0); 1247 } 1248 1249 #undef __FUNCT__ 1250 #define __FUNCT__ "MatSetOption_SeqDense" 1251 PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op,PetscTruth flg) 1252 { 1253 Mat_SeqDense *aij = (Mat_SeqDense*)A->data; 1254 PetscErrorCode ierr; 1255 1256 PetscFunctionBegin; 1257 switch (op) { 1258 case MAT_ROW_ORIENTED: 1259 aij->roworiented = flg; 1260 break; 1261 case MAT_NEW_NONZERO_LOCATIONS: 1262 case MAT_NEW_NONZERO_LOCATION_ERR: 1263 case MAT_NEW_DIAGONALS: 1264 case MAT_IGNORE_OFF_PROC_ENTRIES: 1265 case MAT_USE_HASH_TABLE: 1266 case MAT_SYMMETRIC: 1267 case MAT_STRUCTURALLY_SYMMETRIC: 1268 case MAT_HERMITIAN: 1269 case MAT_SYMMETRY_ETERNAL: 1270 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 1271 break; 1272 default: 1273 SETERRQ1(PETSC_ERR_SUP,"unknown option %d",op); 1274 } 1275 PetscFunctionReturn(0); 1276 } 1277 1278 #undef __FUNCT__ 1279 #define __FUNCT__ "MatZeroEntries_SeqDense" 1280 PetscErrorCode MatZeroEntries_SeqDense(Mat A) 1281 { 1282 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1283 PetscErrorCode ierr; 1284 PetscInt lda=l->lda,m=A->rmap.n,j; 1285 1286 PetscFunctionBegin; 1287 if (lda>m) { 1288 for (j=0; j<A->cmap.n; j++) { 1289 ierr = PetscMemzero(l->v+j*lda,m*sizeof(PetscScalar));CHKERRQ(ierr); 1290 } 1291 } else { 1292 ierr = PetscMemzero(l->v,A->rmap.n*A->cmap.n*sizeof(PetscScalar));CHKERRQ(ierr); 1293 } 1294 PetscFunctionReturn(0); 1295 } 1296 1297 #undef __FUNCT__ 1298 #define __FUNCT__ "MatZeroRows_SeqDense" 1299 PetscErrorCode MatZeroRows_SeqDense(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag) 1300 { 1301 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1302 PetscInt n = A->cmap.n,i,j; 1303 PetscScalar *slot; 1304 1305 PetscFunctionBegin; 1306 for (i=0; i<N; i++) { 1307 slot = l->v + rows[i]; 1308 for (j=0; j<n; j++) { *slot = 0.0; slot += n;} 1309 } 1310 if (diag != 0.0) { 1311 for (i=0; i<N; i++) { 1312 slot = l->v + (n+1)*rows[i]; 1313 *slot = diag; 1314 } 1315 } 1316 PetscFunctionReturn(0); 1317 } 1318 1319 #undef __FUNCT__ 1320 #define __FUNCT__ "MatGetArray_SeqDense" 1321 PetscErrorCode MatGetArray_SeqDense(Mat A,PetscScalar *array[]) 1322 { 1323 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1324 1325 PetscFunctionBegin; 1326 if (mat->lda != A->rmap.n) SETERRQ(PETSC_ERR_SUP,"Cannot get array for Dense matrices with LDA different from number of rows"); 1327 *array = mat->v; 1328 PetscFunctionReturn(0); 1329 } 1330 1331 #undef __FUNCT__ 1332 #define __FUNCT__ "MatRestoreArray_SeqDense" 1333 PetscErrorCode MatRestoreArray_SeqDense(Mat A,PetscScalar *array[]) 1334 { 1335 PetscFunctionBegin; 1336 *array = 0; /* user cannot accidently use the array later */ 1337 PetscFunctionReturn(0); 1338 } 1339 1340 #undef __FUNCT__ 1341 #define __FUNCT__ "MatGetSubMatrix_SeqDense" 1342 static PetscErrorCode MatGetSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B) 1343 { 1344 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1345 PetscErrorCode ierr; 1346 PetscInt i,j,*irow,*icol,nrows,ncols; 1347 PetscScalar *av,*bv,*v = mat->v; 1348 Mat newmat; 1349 1350 PetscFunctionBegin; 1351 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1352 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1353 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1354 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1355 1356 /* Check submatrixcall */ 1357 if (scall == MAT_REUSE_MATRIX) { 1358 PetscInt n_cols,n_rows; 1359 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1360 if (n_rows != nrows || n_cols != ncols) { 1361 /* resize the result result matrix to match number of requested rows/columns */ 1362 ierr = MatSetSizes(*B,nrows,nrows,nrows,nrows);CHKERRQ(ierr); 1363 } 1364 newmat = *B; 1365 } else { 1366 /* Create and fill new matrix */ 1367 ierr = MatCreate(A->comm,&newmat);CHKERRQ(ierr); 1368 ierr = MatSetSizes(newmat,nrows,ncols,nrows,ncols);CHKERRQ(ierr); 1369 ierr = MatSetType(newmat,A->type_name);CHKERRQ(ierr); 1370 ierr = MatSeqDenseSetPreallocation(newmat,PETSC_NULL);CHKERRQ(ierr); 1371 } 1372 1373 /* Now extract the data pointers and do the copy,column at a time */ 1374 bv = ((Mat_SeqDense*)newmat->data)->v; 1375 1376 for (i=0; i<ncols; i++) { 1377 av = v + mat->lda*icol[i]; 1378 for (j=0; j<nrows; j++) { 1379 *bv++ = av[irow[j]]; 1380 } 1381 } 1382 1383 /* Assemble the matrices so that the correct flags are set */ 1384 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1385 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1386 1387 /* Free work space */ 1388 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1389 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1390 *B = newmat; 1391 PetscFunctionReturn(0); 1392 } 1393 1394 #undef __FUNCT__ 1395 #define __FUNCT__ "MatGetSubMatrices_SeqDense" 1396 PetscErrorCode MatGetSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1397 { 1398 PetscErrorCode ierr; 1399 PetscInt i; 1400 1401 PetscFunctionBegin; 1402 if (scall == MAT_INITIAL_MATRIX) { 1403 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1404 } 1405 1406 for (i=0; i<n; i++) { 1407 ierr = MatGetSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1408 } 1409 PetscFunctionReturn(0); 1410 } 1411 1412 #undef __FUNCT__ 1413 #define __FUNCT__ "MatAssemblyBegin_SeqDense" 1414 PetscErrorCode MatAssemblyBegin_SeqDense(Mat mat,MatAssemblyType mode) 1415 { 1416 PetscFunctionBegin; 1417 PetscFunctionReturn(0); 1418 } 1419 1420 #undef __FUNCT__ 1421 #define __FUNCT__ "MatAssemblyEnd_SeqDense" 1422 PetscErrorCode MatAssemblyEnd_SeqDense(Mat mat,MatAssemblyType mode) 1423 { 1424 PetscFunctionBegin; 1425 PetscFunctionReturn(0); 1426 } 1427 1428 #undef __FUNCT__ 1429 #define __FUNCT__ "MatCopy_SeqDense" 1430 PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str) 1431 { 1432 Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense *)B->data; 1433 PetscErrorCode ierr; 1434 PetscInt lda1=a->lda,lda2=b->lda, m=A->rmap.n,n=A->cmap.n, j; 1435 1436 PetscFunctionBegin; 1437 /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */ 1438 if (A->ops->copy != B->ops->copy) { 1439 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1440 PetscFunctionReturn(0); 1441 } 1442 if (m != B->rmap.n || n != B->cmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"size(B) != size(A)"); 1443 if (lda1>m || lda2>m) { 1444 for (j=0; j<n; j++) { 1445 ierr = PetscMemcpy(b->v+j*lda2,a->v+j*lda1,m*sizeof(PetscScalar));CHKERRQ(ierr); 1446 } 1447 } else { 1448 ierr = PetscMemcpy(b->v,a->v,A->rmap.n*A->cmap.n*sizeof(PetscScalar));CHKERRQ(ierr); 1449 } 1450 PetscFunctionReturn(0); 1451 } 1452 1453 #undef __FUNCT__ 1454 #define __FUNCT__ "MatSetUpPreallocation_SeqDense" 1455 PetscErrorCode MatSetUpPreallocation_SeqDense(Mat A) 1456 { 1457 PetscErrorCode ierr; 1458 1459 PetscFunctionBegin; 1460 ierr = MatSeqDenseSetPreallocation(A,0);CHKERRQ(ierr); 1461 PetscFunctionReturn(0); 1462 } 1463 1464 #undef __FUNCT__ 1465 #define __FUNCT__ "MatSetSizes_SeqDense" 1466 PetscErrorCode MatSetSizes_SeqDense(Mat A,PetscInt m,PetscInt n,PetscInt M,PetscInt N) 1467 { 1468 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1469 PetscErrorCode ierr; 1470 PetscFunctionBegin; 1471 /* this will not be called before lda, Mmax, and Nmax have been set */ 1472 m = PetscMax(m,M); 1473 n = PetscMax(n,N); 1474 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); 1475 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); 1476 A->rmap.n = A->rmap.n = m; 1477 A->cmap.n = A->cmap.N = n; 1478 if (a->changelda) a->lda = m; 1479 ierr = PetscMemzero(a->v,m*n*sizeof(PetscScalar));CHKERRQ(ierr); 1480 PetscFunctionReturn(0); 1481 } 1482 1483 /* ----------------------------------------------------------------*/ 1484 #undef __FUNCT__ 1485 #define __FUNCT__ "MatMatMult_SeqDense_SeqDense" 1486 PetscErrorCode MatMatMult_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 1487 { 1488 PetscErrorCode ierr; 1489 1490 PetscFunctionBegin; 1491 if (scall == MAT_INITIAL_MATRIX){ 1492 ierr = MatMatMultSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 1493 } 1494 ierr = MatMatMultNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 1495 PetscFunctionReturn(0); 1496 } 1497 1498 #undef __FUNCT__ 1499 #define __FUNCT__ "MatMatMultSymbolic_SeqDense_SeqDense" 1500 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 1501 { 1502 PetscErrorCode ierr; 1503 PetscInt m=A->rmap.n,n=B->cmap.n; 1504 Mat Cmat; 1505 1506 PetscFunctionBegin; 1507 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); 1508 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 1509 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 1510 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 1511 ierr = MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);CHKERRQ(ierr); 1512 Cmat->assembled = PETSC_TRUE; 1513 *C = Cmat; 1514 PetscFunctionReturn(0); 1515 } 1516 1517 #undef __FUNCT__ 1518 #define __FUNCT__ "MatMatMultNumeric_SeqDense_SeqDense" 1519 PetscErrorCode MatMatMultNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 1520 { 1521 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1522 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1523 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 1524 PetscBLASInt m=(PetscBLASInt)A->rmap.n,n=(PetscBLASInt)B->cmap.n,k=(PetscBLASInt)A->cmap.n; 1525 PetscScalar _DOne=1.0,_DZero=0.0; 1526 1527 PetscFunctionBegin; 1528 BLASgemm_("N","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda); 1529 PetscFunctionReturn(0); 1530 } 1531 1532 #undef __FUNCT__ 1533 #define __FUNCT__ "MatMatMultTranspose_SeqDense_SeqDense" 1534 PetscErrorCode MatMatMultTranspose_SeqDense_SeqDense(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 1535 { 1536 PetscErrorCode ierr; 1537 1538 PetscFunctionBegin; 1539 if (scall == MAT_INITIAL_MATRIX){ 1540 ierr = MatMatMultTransposeSymbolic_SeqDense_SeqDense(A,B,fill,C);CHKERRQ(ierr); 1541 } 1542 ierr = MatMatMultTransposeNumeric_SeqDense_SeqDense(A,B,*C);CHKERRQ(ierr); 1543 PetscFunctionReturn(0); 1544 } 1545 1546 #undef __FUNCT__ 1547 #define __FUNCT__ "MatMatMultTransposeSymbolic_SeqDense_SeqDense" 1548 PetscErrorCode MatMatMultTransposeSymbolic_SeqDense_SeqDense(Mat A,Mat B,PetscReal fill,Mat *C) 1549 { 1550 PetscErrorCode ierr; 1551 PetscInt m=A->cmap.n,n=B->cmap.n; 1552 Mat Cmat; 1553 1554 PetscFunctionBegin; 1555 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); 1556 ierr = MatCreate(PETSC_COMM_SELF,&Cmat);CHKERRQ(ierr); 1557 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 1558 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 1559 ierr = MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);CHKERRQ(ierr); 1560 Cmat->assembled = PETSC_TRUE; 1561 *C = Cmat; 1562 PetscFunctionReturn(0); 1563 } 1564 1565 #undef __FUNCT__ 1566 #define __FUNCT__ "MatMatMultTransposeNumeric_SeqDense_SeqDense" 1567 PetscErrorCode MatMatMultTransposeNumeric_SeqDense_SeqDense(Mat A,Mat B,Mat C) 1568 { 1569 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1570 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1571 Mat_SeqDense *c = (Mat_SeqDense*)C->data; 1572 PetscBLASInt m=(PetscBLASInt)A->cmap.n,n=(PetscBLASInt)B->cmap.n,k=(PetscBLASInt)A->rmap.n; 1573 PetscScalar _DOne=1.0,_DZero=0.0; 1574 1575 PetscFunctionBegin; 1576 BLASgemm_("T","N",&m,&n,&k,&_DOne,a->v,&a->lda,b->v,&b->lda,&_DZero,c->v,&c->lda); 1577 PetscFunctionReturn(0); 1578 } 1579 1580 #undef __FUNCT__ 1581 #define __FUNCT__ "MatGetRowMax_SeqDense" 1582 PetscErrorCode MatGetRowMax_SeqDense(Mat A,Vec v,PetscInt idx[]) 1583 { 1584 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1585 PetscErrorCode ierr; 1586 PetscInt i,j,m = A->rmap.n,n = A->cmap.n,p; 1587 PetscScalar *x; 1588 MatScalar *aa = a->v; 1589 1590 PetscFunctionBegin; 1591 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1592 1593 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1594 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1595 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1596 if (p != A->rmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1597 for (i=0; i<m; i++) { 1598 x[i] = aa[i]; if (idx) idx[i] = 0; 1599 for (j=1; j<n; j++){ 1600 if (PetscRealPart(x[i]) < PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;} 1601 } 1602 } 1603 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1604 PetscFunctionReturn(0); 1605 } 1606 1607 #undef __FUNCT__ 1608 #define __FUNCT__ "MatGetRowMaxAbs_SeqDense" 1609 PetscErrorCode MatGetRowMaxAbs_SeqDense(Mat A,Vec v,PetscInt idx[]) 1610 { 1611 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1612 PetscErrorCode ierr; 1613 PetscInt i,j,m = A->rmap.n,n = A->cmap.n,p; 1614 PetscScalar *x; 1615 PetscReal atmp; 1616 MatScalar *aa = a->v; 1617 1618 PetscFunctionBegin; 1619 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1620 1621 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1622 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1623 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1624 if (p != A->rmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1625 for (i=0; i<m; i++) { 1626 x[i] = PetscAbsScalar(aa[i]); if (idx) idx[i] = 0; 1627 for (j=1; j<n; j++){ 1628 atmp = PetscAbsScalar(aa[i+m*j]); 1629 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = j;} 1630 } 1631 } 1632 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1633 PetscFunctionReturn(0); 1634 } 1635 1636 #undef __FUNCT__ 1637 #define __FUNCT__ "MatGetRowMin_SeqDense" 1638 PetscErrorCode MatGetRowMin_SeqDense(Mat A,Vec v,PetscInt idx[]) 1639 { 1640 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1641 PetscErrorCode ierr; 1642 PetscInt i,j,m = A->rmap.n,n = A->cmap.n,p; 1643 PetscScalar *x; 1644 MatScalar *aa = a->v; 1645 1646 PetscFunctionBegin; 1647 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1648 1649 ierr = VecSet(v,0.0);CHKERRQ(ierr); 1650 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1651 ierr = VecGetLocalSize(v,&p);CHKERRQ(ierr); 1652 if (p != A->rmap.n) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1653 for (i=0; i<m; i++) { 1654 x[i] = aa[i]; if (idx) idx[i] = 0; 1655 for (j=1; j<n; j++){ 1656 if (PetscRealPart(x[i]) > PetscRealPart(aa[i+m*j])) {x[i] = aa[i + m*j]; if (idx) idx[i] = j;} 1657 } 1658 } 1659 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1660 PetscFunctionReturn(0); 1661 } 1662 1663 #undef __FUNCT__ 1664 #define __FUNCT__ "MatGetColumnVector_SeqDense" 1665 PetscErrorCode MatGetColumnVector_SeqDense(Mat A,Vec v,PetscInt col) 1666 { 1667 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 1668 PetscErrorCode ierr; 1669 PetscScalar *x; 1670 1671 PetscFunctionBegin; 1672 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 1673 1674 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1675 ierr = PetscMemcpy(x,a->v+col*a->lda,A->rmap.n*sizeof(PetscScalar));CHKERRQ(ierr); 1676 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1677 PetscFunctionReturn(0); 1678 } 1679 1680 /* -------------------------------------------------------------------*/ 1681 static struct _MatOps MatOps_Values = {MatSetValues_SeqDense, 1682 MatGetRow_SeqDense, 1683 MatRestoreRow_SeqDense, 1684 MatMult_SeqDense, 1685 /* 4*/ MatMultAdd_SeqDense, 1686 MatMultTranspose_SeqDense, 1687 MatMultTransposeAdd_SeqDense, 1688 MatSolve_SeqDense, 1689 MatSolveAdd_SeqDense, 1690 MatSolveTranspose_SeqDense, 1691 /*10*/ MatSolveTransposeAdd_SeqDense, 1692 MatLUFactor_SeqDense, 1693 MatCholeskyFactor_SeqDense, 1694 MatRelax_SeqDense, 1695 MatTranspose_SeqDense, 1696 /*15*/ MatGetInfo_SeqDense, 1697 MatEqual_SeqDense, 1698 MatGetDiagonal_SeqDense, 1699 MatDiagonalScale_SeqDense, 1700 MatNorm_SeqDense, 1701 /*20*/ MatAssemblyBegin_SeqDense, 1702 MatAssemblyEnd_SeqDense, 1703 0, 1704 MatSetOption_SeqDense, 1705 MatZeroEntries_SeqDense, 1706 /*25*/ MatZeroRows_SeqDense, 1707 MatLUFactorSymbolic_SeqDense, 1708 MatLUFactorNumeric_SeqDense, 1709 MatCholeskyFactorSymbolic_SeqDense, 1710 MatCholeskyFactorNumeric_SeqDense, 1711 /*30*/ MatSetUpPreallocation_SeqDense, 1712 0, 1713 0, 1714 MatGetArray_SeqDense, 1715 MatRestoreArray_SeqDense, 1716 /*35*/ MatDuplicate_SeqDense, 1717 0, 1718 0, 1719 0, 1720 0, 1721 /*40*/ MatAXPY_SeqDense, 1722 MatGetSubMatrices_SeqDense, 1723 0, 1724 MatGetValues_SeqDense, 1725 MatCopy_SeqDense, 1726 /*45*/ MatGetRowMax_SeqDense, 1727 MatScale_SeqDense, 1728 0, 1729 0, 1730 0, 1731 /*50*/ 0, 1732 0, 1733 0, 1734 0, 1735 0, 1736 /*55*/ 0, 1737 0, 1738 0, 1739 0, 1740 0, 1741 /*60*/ 0, 1742 MatDestroy_SeqDense, 1743 MatView_SeqDense, 1744 0, 1745 0, 1746 /*65*/ 0, 1747 0, 1748 0, 1749 0, 1750 0, 1751 /*70*/ MatGetRowMaxAbs_SeqDense, 1752 0, 1753 0, 1754 0, 1755 0, 1756 /*75*/ 0, 1757 0, 1758 0, 1759 0, 1760 0, 1761 /*80*/ 0, 1762 0, 1763 0, 1764 0, 1765 /*84*/ MatLoad_SeqDense, 1766 0, 1767 MatIsHermitian_SeqDense, 1768 0, 1769 0, 1770 0, 1771 /*90*/ MatMatMult_SeqDense_SeqDense, 1772 MatMatMultSymbolic_SeqDense_SeqDense, 1773 MatMatMultNumeric_SeqDense_SeqDense, 1774 0, 1775 0, 1776 /*95*/ 0, 1777 MatMatMultTranspose_SeqDense_SeqDense, 1778 MatMatMultTransposeSymbolic_SeqDense_SeqDense, 1779 MatMatMultTransposeNumeric_SeqDense_SeqDense, 1780 0, 1781 /*100*/0, 1782 0, 1783 0, 1784 0, 1785 MatSetSizes_SeqDense, 1786 0, 1787 0, 1788 0, 1789 0, 1790 0, 1791 /*110*/0, 1792 0, 1793 MatGetRowMin_SeqDense, 1794 MatGetColumnVector_SeqDense 1795 }; 1796 1797 #undef __FUNCT__ 1798 #define __FUNCT__ "MatCreateSeqDense" 1799 /*@C 1800 MatCreateSeqDense - Creates a sequential dense matrix that 1801 is stored in column major order (the usual Fortran 77 manner). Many 1802 of the matrix operations use the BLAS and LAPACK routines. 1803 1804 Collective on MPI_Comm 1805 1806 Input Parameters: 1807 + comm - MPI communicator, set to PETSC_COMM_SELF 1808 . m - number of rows 1809 . n - number of columns 1810 - data - optional location of matrix data. Set data=PETSC_NULL for PETSc 1811 to control all matrix memory allocation. 1812 1813 Output Parameter: 1814 . A - the matrix 1815 1816 Notes: 1817 The data input variable is intended primarily for Fortran programmers 1818 who wish to allocate their own matrix memory space. Most users should 1819 set data=PETSC_NULL. 1820 1821 Level: intermediate 1822 1823 .keywords: dense, matrix, LAPACK, BLAS 1824 1825 .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues() 1826 @*/ 1827 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A) 1828 { 1829 PetscErrorCode ierr; 1830 1831 PetscFunctionBegin; 1832 ierr = MatCreate(comm,A);CHKERRQ(ierr); 1833 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 1834 ierr = MatSetType(*A,MATSEQDENSE);CHKERRQ(ierr); 1835 ierr = MatSeqDenseSetPreallocation(*A,data);CHKERRQ(ierr); 1836 PetscFunctionReturn(0); 1837 } 1838 1839 #undef __FUNCT__ 1840 #define __FUNCT__ "MatSeqDenseSetPreallocation" 1841 /*@C 1842 MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements 1843 1844 Collective on MPI_Comm 1845 1846 Input Parameters: 1847 + A - the matrix 1848 - data - the array (or PETSC_NULL) 1849 1850 Notes: 1851 The data input variable is intended primarily for Fortran programmers 1852 who wish to allocate their own matrix memory space. Most users should 1853 need not call this routine. 1854 1855 Level: intermediate 1856 1857 .keywords: dense, matrix, LAPACK, BLAS 1858 1859 .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues() 1860 @*/ 1861 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqDenseSetPreallocation(Mat B,PetscScalar data[]) 1862 { 1863 PetscErrorCode ierr,(*f)(Mat,PetscScalar[]); 1864 1865 PetscFunctionBegin; 1866 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 1867 if (f) { 1868 ierr = (*f)(B,data);CHKERRQ(ierr); 1869 } 1870 PetscFunctionReturn(0); 1871 } 1872 1873 EXTERN_C_BEGIN 1874 #undef __FUNCT__ 1875 #define __FUNCT__ "MatSeqDenseSetPreallocation_SeqDense" 1876 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data) 1877 { 1878 Mat_SeqDense *b; 1879 PetscErrorCode ierr; 1880 1881 PetscFunctionBegin; 1882 B->preallocated = PETSC_TRUE; 1883 b = (Mat_SeqDense*)B->data; 1884 if (!data) { 1885 ierr = PetscMalloc((b->lda*b->Nmax+1)*sizeof(PetscScalar),&b->v);CHKERRQ(ierr); 1886 ierr = PetscMemzero(b->v,b->lda*b->Nmax*sizeof(PetscScalar));CHKERRQ(ierr); 1887 b->user_alloc = PETSC_FALSE; 1888 ierr = PetscLogObjectMemory(B,b->lda*b->Nmax*sizeof(PetscScalar));CHKERRQ(ierr); 1889 } else { /* user-allocated storage */ 1890 b->v = data; 1891 b->user_alloc = PETSC_TRUE; 1892 } 1893 PetscFunctionReturn(0); 1894 } 1895 EXTERN_C_END 1896 1897 #undef __FUNCT__ 1898 #define __FUNCT__ "MatSeqDenseSetLDA" 1899 /*@C 1900 MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array 1901 1902 Input parameter: 1903 + A - the matrix 1904 - lda - the leading dimension 1905 1906 Notes: 1907 This routine is to be used in conjunction with MatSeqDenseSetPreallocation; 1908 it asserts that the preallocation has a leading dimension (the LDA parameter 1909 of Blas and Lapack fame) larger than M, the first dimension of the matrix. 1910 1911 Level: intermediate 1912 1913 .keywords: dense, matrix, LAPACK, BLAS 1914 1915 .seealso: MatCreate(), MatCreateSeqDense(), MatSeqDenseSetPreallocation(), MatSetMaximumSize() 1916 @*/ 1917 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqDenseSetLDA(Mat B,PetscInt lda) 1918 { 1919 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1920 1921 PetscFunctionBegin; 1922 if (lda < B->rmap.n) SETERRQ2(PETSC_ERR_ARG_SIZ,"LDA %D must be at least matrix dimension %D",lda,B->rmap.n); 1923 b->lda = lda; 1924 b->changelda = PETSC_FALSE; 1925 b->Mmax = PetscMax(b->Mmax,lda); 1926 PetscFunctionReturn(0); 1927 } 1928 1929 /*MC 1930 MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices. 1931 1932 Options Database Keys: 1933 . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions() 1934 1935 Level: beginner 1936 1937 .seealso: MatCreateSeqDense() 1938 1939 M*/ 1940 1941 EXTERN_C_BEGIN 1942 #undef __FUNCT__ 1943 #define __FUNCT__ "MatCreate_SeqDense" 1944 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_SeqDense(Mat B) 1945 { 1946 Mat_SeqDense *b; 1947 PetscErrorCode ierr; 1948 PetscMPIInt size; 1949 1950 PetscFunctionBegin; 1951 ierr = MPI_Comm_size(B->comm,&size);CHKERRQ(ierr); 1952 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 1953 1954 B->rmap.bs = B->cmap.bs = 1; 1955 ierr = PetscMapSetUp(&B->rmap);CHKERRQ(ierr); 1956 ierr = PetscMapSetUp(&B->cmap);CHKERRQ(ierr); 1957 1958 ierr = PetscNewLog(B,Mat_SeqDense,&b);CHKERRQ(ierr); 1959 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1960 B->factor = 0; 1961 B->mapping = 0; 1962 B->data = (void*)b; 1963 1964 1965 b->pivots = 0; 1966 b->roworiented = PETSC_TRUE; 1967 b->v = 0; 1968 b->lda = B->rmap.n; 1969 b->changelda = PETSC_FALSE; 1970 b->Mmax = B->rmap.n; 1971 b->Nmax = B->cmap.n; 1972 1973 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqDenseSetPreallocation_C", 1974 "MatSeqDenseSetPreallocation_SeqDense", 1975 MatSeqDenseSetPreallocation_SeqDense);CHKERRQ(ierr); 1976 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMult_seqaij_seqdense_C", 1977 "MatMatMult_SeqAIJ_SeqDense", 1978 MatMatMult_SeqAIJ_SeqDense);CHKERRQ(ierr); 1979 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultSymbolic_seqaij_seqdense_C", 1980 "MatMatMultSymbolic_SeqAIJ_SeqDense", 1981 MatMatMultSymbolic_SeqAIJ_SeqDense);CHKERRQ(ierr); 1982 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultNumeric_seqaij_seqdense_C", 1983 "MatMatMultNumeric_SeqAIJ_SeqDense", 1984 MatMatMultNumeric_SeqAIJ_SeqDense);CHKERRQ(ierr); 1985 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQDENSE);CHKERRQ(ierr); 1986 PetscFunctionReturn(0); 1987 } 1988 1989 1990 EXTERN_C_END 1991