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