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