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