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