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