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