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 PetscInt 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,PETSC_TRUE);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,PETSC_FALSE);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,PETSC_TRUE);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,PETSC_FALSE);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,PETSC_FALSE);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 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatSeqDenseSetPreallocation_C","",PETSC_NULL);CHKERRQ(ierr); 1031 PetscFunctionReturn(0); 1032 } 1033 1034 #undef __FUNCT__ 1035 #define __FUNCT__ "MatTranspose_SeqDense" 1036 PetscErrorCode MatTranspose_SeqDense(Mat A,Mat *matout) 1037 { 1038 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1039 PetscErrorCode ierr; 1040 int k,j,m,n,M; 1041 PetscScalar *v,tmp; 1042 1043 PetscFunctionBegin; 1044 v = mat->v; m = A->m; M = mat->lda; n = A->n; 1045 if (!matout) { /* in place transpose */ 1046 if (m != n) { 1047 SETERRQ(PETSC_ERR_SUP,"Can not transpose non-square matrix in place"); 1048 } else { 1049 for (j=0; j<m; j++) { 1050 for (k=0; k<j; k++) { 1051 tmp = v[j + k*M]; 1052 v[j + k*M] = v[k + j*M]; 1053 v[k + j*M] = tmp; 1054 } 1055 } 1056 } 1057 } else { /* out-of-place transpose */ 1058 Mat tmat; 1059 Mat_SeqDense *tmatd; 1060 PetscScalar *v2; 1061 1062 ierr = MatCreate(A->comm,A->n,A->m,A->n,A->m,&tmat);CHKERRQ(ierr); 1063 ierr = MatSetType(tmat,A->type_name);CHKERRQ(ierr); 1064 ierr = MatSeqDenseSetPreallocation(tmat,PETSC_NULL);CHKERRQ(ierr); 1065 tmatd = (Mat_SeqDense*)tmat->data; 1066 v = mat->v; v2 = tmatd->v; 1067 for (j=0; j<n; j++) { 1068 for (k=0; k<m; k++) v2[j + k*n] = v[k + j*M]; 1069 } 1070 ierr = MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1071 ierr = MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1072 *matout = tmat; 1073 } 1074 PetscFunctionReturn(0); 1075 } 1076 1077 #undef __FUNCT__ 1078 #define __FUNCT__ "MatEqual_SeqDense" 1079 PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscTruth *flg) 1080 { 1081 Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data; 1082 Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data; 1083 int i,j; 1084 PetscScalar *v1 = mat1->v,*v2 = mat2->v; 1085 1086 PetscFunctionBegin; 1087 if (A1->m != A2->m) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1088 if (A1->n != A2->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1089 for (i=0; i<A1->m; i++) { 1090 v1 = mat1->v+i; v2 = mat2->v+i; 1091 for (j=0; j<A1->n; j++) { 1092 if (*v1 != *v2) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1093 v1 += mat1->lda; v2 += mat2->lda; 1094 } 1095 } 1096 *flg = PETSC_TRUE; 1097 PetscFunctionReturn(0); 1098 } 1099 1100 #undef __FUNCT__ 1101 #define __FUNCT__ "MatGetDiagonal_SeqDense" 1102 PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v) 1103 { 1104 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1105 PetscErrorCode ierr; 1106 int i,n,len; 1107 PetscScalar *x,zero = 0.0; 1108 1109 PetscFunctionBegin; 1110 ierr = VecSet(&zero,v);CHKERRQ(ierr); 1111 ierr = VecGetSize(v,&n);CHKERRQ(ierr); 1112 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1113 len = PetscMin(A->m,A->n); 1114 if (n != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec"); 1115 for (i=0; i<len; i++) { 1116 x[i] = mat->v[i*mat->lda + i]; 1117 } 1118 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1119 PetscFunctionReturn(0); 1120 } 1121 1122 #undef __FUNCT__ 1123 #define __FUNCT__ "MatDiagonalScale_SeqDense" 1124 PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr) 1125 { 1126 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1127 PetscScalar *l,*r,x,*v; 1128 PetscErrorCode ierr; 1129 int i,j,m = A->m,n = A->n; 1130 1131 PetscFunctionBegin; 1132 if (ll) { 1133 ierr = VecGetSize(ll,&m);CHKERRQ(ierr); 1134 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1135 if (m != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size"); 1136 for (i=0; i<m; i++) { 1137 x = l[i]; 1138 v = mat->v + i; 1139 for (j=0; j<n; j++) { (*v) *= x; v+= m;} 1140 } 1141 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1142 PetscLogFlops(n*m); 1143 } 1144 if (rr) { 1145 ierr = VecGetSize(rr,&n);CHKERRQ(ierr); 1146 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1147 if (n != A->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size"); 1148 for (i=0; i<n; i++) { 1149 x = r[i]; 1150 v = mat->v + i*m; 1151 for (j=0; j<m; j++) { (*v++) *= x;} 1152 } 1153 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1154 PetscLogFlops(n*m); 1155 } 1156 PetscFunctionReturn(0); 1157 } 1158 1159 #undef __FUNCT__ 1160 #define __FUNCT__ "MatNorm_SeqDense" 1161 PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm) 1162 { 1163 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1164 PetscScalar *v = mat->v; 1165 PetscReal sum = 0.0; 1166 int lda=mat->lda,m=A->m,i,j; 1167 1168 PetscFunctionBegin; 1169 if (type == NORM_FROBENIUS) { 1170 if (lda>m) { 1171 for (j=0; j<A->n; j++) { 1172 v = mat->v+j*lda; 1173 for (i=0; i<m; i++) { 1174 #if defined(PETSC_USE_COMPLEX) 1175 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1176 #else 1177 sum += (*v)*(*v); v++; 1178 #endif 1179 } 1180 } 1181 } else { 1182 for (i=0; i<A->n*A->m; i++) { 1183 #if defined(PETSC_USE_COMPLEX) 1184 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1185 #else 1186 sum += (*v)*(*v); v++; 1187 #endif 1188 } 1189 } 1190 *nrm = sqrt(sum); 1191 PetscLogFlops(2*A->n*A->m); 1192 } else if (type == NORM_1) { 1193 *nrm = 0.0; 1194 for (j=0; j<A->n; j++) { 1195 v = mat->v + j*mat->lda; 1196 sum = 0.0; 1197 for (i=0; i<A->m; i++) { 1198 sum += PetscAbsScalar(*v); v++; 1199 } 1200 if (sum > *nrm) *nrm = sum; 1201 } 1202 PetscLogFlops(A->n*A->m); 1203 } else if (type == NORM_INFINITY) { 1204 *nrm = 0.0; 1205 for (j=0; j<A->m; j++) { 1206 v = mat->v + j; 1207 sum = 0.0; 1208 for (i=0; i<A->n; i++) { 1209 sum += PetscAbsScalar(*v); v += mat->lda; 1210 } 1211 if (sum > *nrm) *nrm = sum; 1212 } 1213 PetscLogFlops(A->n*A->m); 1214 } else { 1215 SETERRQ(PETSC_ERR_SUP,"No two norm"); 1216 } 1217 PetscFunctionReturn(0); 1218 } 1219 1220 #undef __FUNCT__ 1221 #define __FUNCT__ "MatSetOption_SeqDense" 1222 PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op) 1223 { 1224 Mat_SeqDense *aij = (Mat_SeqDense*)A->data; 1225 1226 PetscFunctionBegin; 1227 switch (op) { 1228 case MAT_ROW_ORIENTED: 1229 aij->roworiented = PETSC_TRUE; 1230 break; 1231 case MAT_COLUMN_ORIENTED: 1232 aij->roworiented = PETSC_FALSE; 1233 break; 1234 case MAT_ROWS_SORTED: 1235 case MAT_ROWS_UNSORTED: 1236 case MAT_COLUMNS_SORTED: 1237 case MAT_COLUMNS_UNSORTED: 1238 case MAT_NO_NEW_NONZERO_LOCATIONS: 1239 case MAT_YES_NEW_NONZERO_LOCATIONS: 1240 case MAT_NEW_NONZERO_LOCATION_ERR: 1241 case MAT_NO_NEW_DIAGONALS: 1242 case MAT_YES_NEW_DIAGONALS: 1243 case MAT_IGNORE_OFF_PROC_ENTRIES: 1244 case MAT_USE_HASH_TABLE: 1245 PetscLogInfo(A,"MatSetOption_SeqDense:Option ignored\n"); 1246 break; 1247 case MAT_SYMMETRIC: 1248 case MAT_STRUCTURALLY_SYMMETRIC: 1249 case MAT_NOT_SYMMETRIC: 1250 case MAT_NOT_STRUCTURALLY_SYMMETRIC: 1251 case MAT_HERMITIAN: 1252 case MAT_NOT_HERMITIAN: 1253 case MAT_SYMMETRY_ETERNAL: 1254 case MAT_NOT_SYMMETRY_ETERNAL: 1255 break; 1256 default: 1257 SETERRQ(PETSC_ERR_SUP,"unknown option"); 1258 } 1259 PetscFunctionReturn(0); 1260 } 1261 1262 #undef __FUNCT__ 1263 #define __FUNCT__ "MatZeroEntries_SeqDense" 1264 PetscErrorCode MatZeroEntries_SeqDense(Mat A) 1265 { 1266 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1267 PetscErrorCode ierr; 1268 int lda=l->lda,m=A->m,j; 1269 1270 PetscFunctionBegin; 1271 if (lda>m) { 1272 for (j=0; j<A->n; j++) { 1273 ierr = PetscMemzero(l->v+j*lda,m*sizeof(PetscScalar));CHKERRQ(ierr); 1274 } 1275 } else { 1276 ierr = PetscMemzero(l->v,A->m*A->n*sizeof(PetscScalar));CHKERRQ(ierr); 1277 } 1278 PetscFunctionReturn(0); 1279 } 1280 1281 #undef __FUNCT__ 1282 #define __FUNCT__ "MatGetBlockSize_SeqDense" 1283 PetscErrorCode MatGetBlockSize_SeqDense(Mat A,int *bs) 1284 { 1285 PetscFunctionBegin; 1286 *bs = 1; 1287 PetscFunctionReturn(0); 1288 } 1289 1290 #undef __FUNCT__ 1291 #define __FUNCT__ "MatZeroRows_SeqDense" 1292 PetscErrorCode MatZeroRows_SeqDense(Mat A,IS is,const PetscScalar *diag) 1293 { 1294 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1295 PetscErrorCode ierr; 1296 int n = A->n,i,j,N,*rows; 1297 PetscScalar *slot; 1298 1299 PetscFunctionBegin; 1300 ierr = ISGetLocalSize(is,&N);CHKERRQ(ierr); 1301 ierr = ISGetIndices(is,&rows);CHKERRQ(ierr); 1302 for (i=0; i<N; i++) { 1303 slot = l->v + rows[i]; 1304 for (j=0; j<n; j++) { *slot = 0.0; slot += n;} 1305 } 1306 if (diag) { 1307 for (i=0; i<N; i++) { 1308 slot = l->v + (n+1)*rows[i]; 1309 *slot = *diag; 1310 } 1311 } 1312 ISRestoreIndices(is,&rows); 1313 PetscFunctionReturn(0); 1314 } 1315 1316 #undef __FUNCT__ 1317 #define __FUNCT__ "MatGetArray_SeqDense" 1318 PetscErrorCode MatGetArray_SeqDense(Mat A,PetscScalar *array[]) 1319 { 1320 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1321 1322 PetscFunctionBegin; 1323 *array = mat->v; 1324 PetscFunctionReturn(0); 1325 } 1326 1327 #undef __FUNCT__ 1328 #define __FUNCT__ "MatRestoreArray_SeqDense" 1329 PetscErrorCode MatRestoreArray_SeqDense(Mat A,PetscScalar *array[]) 1330 { 1331 PetscFunctionBegin; 1332 *array = 0; /* user cannot accidently use the array later */ 1333 PetscFunctionReturn(0); 1334 } 1335 1336 #undef __FUNCT__ 1337 #define __FUNCT__ "MatGetSubMatrix_SeqDense" 1338 static PetscErrorCode MatGetSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,int cs,MatReuse scall,Mat *B) 1339 { 1340 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1341 PetscErrorCode ierr; 1342 int i,j,m = A->m,*irow,*icol,nrows,ncols; 1343 PetscScalar *av,*bv,*v = mat->v; 1344 Mat newmat; 1345 1346 PetscFunctionBegin; 1347 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1348 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1349 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1350 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1351 1352 /* Check submatrixcall */ 1353 if (scall == MAT_REUSE_MATRIX) { 1354 int n_cols,n_rows; 1355 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1356 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 1357 newmat = *B; 1358 } else { 1359 /* Create and fill new matrix */ 1360 ierr = MatCreate(A->comm,nrows,ncols,nrows,ncols,&newmat);CHKERRQ(ierr); 1361 ierr = MatSetType(newmat,A->type_name);CHKERRQ(ierr); 1362 ierr = MatSeqDenseSetPreallocation(newmat,PETSC_NULL);CHKERRQ(ierr); 1363 } 1364 1365 /* Now extract the data pointers and do the copy,column at a time */ 1366 bv = ((Mat_SeqDense*)newmat->data)->v; 1367 1368 for (i=0; i<ncols; i++) { 1369 av = v + m*icol[i]; 1370 for (j=0; j<nrows; j++) { 1371 *bv++ = av[irow[j]]; 1372 } 1373 } 1374 1375 /* Assemble the matrices so that the correct flags are set */ 1376 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1377 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1378 1379 /* Free work space */ 1380 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1381 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1382 *B = newmat; 1383 PetscFunctionReturn(0); 1384 } 1385 1386 #undef __FUNCT__ 1387 #define __FUNCT__ "MatGetSubMatrices_SeqDense" 1388 PetscErrorCode MatGetSubMatrices_SeqDense(Mat A,int n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1389 { 1390 PetscErrorCode ierr; 1391 int i; 1392 1393 PetscFunctionBegin; 1394 if (scall == MAT_INITIAL_MATRIX) { 1395 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1396 } 1397 1398 for (i=0; i<n; i++) { 1399 ierr = MatGetSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1400 } 1401 PetscFunctionReturn(0); 1402 } 1403 1404 #undef __FUNCT__ 1405 #define __FUNCT__ "MatCopy_SeqDense" 1406 PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str) 1407 { 1408 Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense *)B->data; 1409 PetscErrorCode ierr; 1410 int lda1=a->lda,lda2=b->lda, m=A->m,n=A->n, j; 1411 1412 PetscFunctionBegin; 1413 /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */ 1414 if (A->ops->copy != B->ops->copy) { 1415 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1416 PetscFunctionReturn(0); 1417 } 1418 if (m != B->m || n != B->n) SETERRQ(PETSC_ERR_ARG_SIZ,"size(B) != size(A)"); 1419 if (lda1>m || lda2>m) { 1420 for (j=0; j<n; j++) { 1421 ierr = PetscMemcpy(b->v+j*lda2,a->v+j*lda1,m*sizeof(PetscScalar));CHKERRQ(ierr); 1422 } 1423 } else { 1424 ierr = PetscMemcpy(b->v,a->v,A->m*A->n*sizeof(PetscScalar));CHKERRQ(ierr); 1425 } 1426 PetscFunctionReturn(0); 1427 } 1428 1429 #undef __FUNCT__ 1430 #define __FUNCT__ "MatSetUpPreallocation_SeqDense" 1431 PetscErrorCode MatSetUpPreallocation_SeqDense(Mat A) 1432 { 1433 PetscErrorCode ierr; 1434 1435 PetscFunctionBegin; 1436 ierr = MatSeqDenseSetPreallocation(A,0);CHKERRQ(ierr); 1437 PetscFunctionReturn(0); 1438 } 1439 1440 /* -------------------------------------------------------------------*/ 1441 static struct _MatOps MatOps_Values = {MatSetValues_SeqDense, 1442 MatGetRow_SeqDense, 1443 MatRestoreRow_SeqDense, 1444 MatMult_SeqDense, 1445 /* 4*/ MatMultAdd_SeqDense, 1446 MatMultTranspose_SeqDense, 1447 MatMultTransposeAdd_SeqDense, 1448 MatSolve_SeqDense, 1449 MatSolveAdd_SeqDense, 1450 MatSolveTranspose_SeqDense, 1451 /*10*/ MatSolveTransposeAdd_SeqDense, 1452 MatLUFactor_SeqDense, 1453 MatCholeskyFactor_SeqDense, 1454 MatRelax_SeqDense, 1455 MatTranspose_SeqDense, 1456 /*15*/ MatGetInfo_SeqDense, 1457 MatEqual_SeqDense, 1458 MatGetDiagonal_SeqDense, 1459 MatDiagonalScale_SeqDense, 1460 MatNorm_SeqDense, 1461 /*20*/ 0, 1462 0, 1463 0, 1464 MatSetOption_SeqDense, 1465 MatZeroEntries_SeqDense, 1466 /*25*/ MatZeroRows_SeqDense, 1467 MatLUFactorSymbolic_SeqDense, 1468 MatLUFactorNumeric_SeqDense, 1469 MatCholeskyFactorSymbolic_SeqDense, 1470 MatCholeskyFactorNumeric_SeqDense, 1471 /*30*/ MatSetUpPreallocation_SeqDense, 1472 0, 1473 0, 1474 MatGetArray_SeqDense, 1475 MatRestoreArray_SeqDense, 1476 /*35*/ MatDuplicate_SeqDense, 1477 0, 1478 0, 1479 0, 1480 0, 1481 /*40*/ MatAXPY_SeqDense, 1482 MatGetSubMatrices_SeqDense, 1483 0, 1484 MatGetValues_SeqDense, 1485 MatCopy_SeqDense, 1486 /*45*/ 0, 1487 MatScale_SeqDense, 1488 0, 1489 0, 1490 0, 1491 /*50*/ MatGetBlockSize_SeqDense, 1492 0, 1493 0, 1494 0, 1495 0, 1496 /*55*/ 0, 1497 0, 1498 0, 1499 0, 1500 0, 1501 /*60*/ 0, 1502 MatDestroy_SeqDense, 1503 MatView_SeqDense, 1504 MatGetPetscMaps_Petsc, 1505 0, 1506 /*65*/ 0, 1507 0, 1508 0, 1509 0, 1510 0, 1511 /*70*/ 0, 1512 0, 1513 0, 1514 0, 1515 0, 1516 /*75*/ 0, 1517 0, 1518 0, 1519 0, 1520 0, 1521 /*80*/ 0, 1522 0, 1523 0, 1524 0, 1525 /*84*/ MatLoad_SeqDense, 1526 0, 1527 0, 1528 0, 1529 0, 1530 0, 1531 /*90*/ 0, 1532 0, 1533 0, 1534 0, 1535 0, 1536 /*95*/ 0, 1537 0, 1538 0, 1539 0}; 1540 1541 #undef __FUNCT__ 1542 #define __FUNCT__ "MatCreateSeqDense" 1543 /*@C 1544 MatCreateSeqDense - Creates a sequential dense matrix that 1545 is stored in column major order (the usual Fortran 77 manner). Many 1546 of the matrix operations use the BLAS and LAPACK routines. 1547 1548 Collective on MPI_Comm 1549 1550 Input Parameters: 1551 + comm - MPI communicator, set to PETSC_COMM_SELF 1552 . m - number of rows 1553 . n - number of columns 1554 - data - optional location of matrix data. Set data=PETSC_NULL for PETSc 1555 to control all matrix memory allocation. 1556 1557 Output Parameter: 1558 . A - the matrix 1559 1560 Notes: 1561 The data input variable is intended primarily for Fortran programmers 1562 who wish to allocate their own matrix memory space. Most users should 1563 set data=PETSC_NULL. 1564 1565 Level: intermediate 1566 1567 .keywords: dense, matrix, LAPACK, BLAS 1568 1569 .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues() 1570 @*/ 1571 PetscErrorCode MatCreateSeqDense(MPI_Comm comm,int m,int n,PetscScalar *data,Mat *A) 1572 { 1573 PetscErrorCode ierr; 1574 1575 PetscFunctionBegin; 1576 ierr = MatCreate(comm,m,n,m,n,A);CHKERRQ(ierr); 1577 ierr = MatSetType(*A,MATSEQDENSE);CHKERRQ(ierr); 1578 ierr = MatSeqDenseSetPreallocation(*A,data);CHKERRQ(ierr); 1579 PetscFunctionReturn(0); 1580 } 1581 1582 #undef __FUNCT__ 1583 #define __FUNCT__ "MatSeqDensePreallocation" 1584 /*@C 1585 MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements 1586 1587 Collective on MPI_Comm 1588 1589 Input Parameters: 1590 + A - the matrix 1591 - data - the array (or PETSC_NULL) 1592 1593 Notes: 1594 The data input variable is intended primarily for Fortran programmers 1595 who wish to allocate their own matrix memory space. Most users should 1596 set data=PETSC_NULL. 1597 1598 Level: intermediate 1599 1600 .keywords: dense, matrix, LAPACK, BLAS 1601 1602 .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues() 1603 @*/ 1604 PetscErrorCode MatSeqDenseSetPreallocation(Mat B,PetscScalar data[]) 1605 { 1606 PetscErrorCode ierr,(*f)(Mat,PetscScalar[]); 1607 1608 PetscFunctionBegin; 1609 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 1610 if (f) { 1611 ierr = (*f)(B,data);CHKERRQ(ierr); 1612 } 1613 PetscFunctionReturn(0); 1614 } 1615 1616 EXTERN_C_BEGIN 1617 #undef __FUNCT__ 1618 #define __FUNCT__ "MatSeqDensePreallocation_SeqDense" 1619 PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data) 1620 { 1621 Mat_SeqDense *b; 1622 PetscErrorCode ierr; 1623 1624 PetscFunctionBegin; 1625 B->preallocated = PETSC_TRUE; 1626 b = (Mat_SeqDense*)B->data; 1627 if (!data) { 1628 ierr = PetscMalloc((B->m*B->n+1)*sizeof(PetscScalar),&b->v);CHKERRQ(ierr); 1629 ierr = PetscMemzero(b->v,B->m*B->n*sizeof(PetscScalar));CHKERRQ(ierr); 1630 b->user_alloc = PETSC_FALSE; 1631 PetscLogObjectMemory(B,B->n*B->m*sizeof(PetscScalar)); 1632 } else { /* user-allocated storage */ 1633 b->v = data; 1634 b->user_alloc = PETSC_TRUE; 1635 } 1636 PetscFunctionReturn(0); 1637 } 1638 EXTERN_C_END 1639 1640 #undef __FUNCT__ 1641 #define __FUNCT__ "MatSeqDenseSetLDA" 1642 /*@C 1643 MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array 1644 1645 Input parameter: 1646 + A - the matrix 1647 - lda - the leading dimension 1648 1649 Notes: 1650 This routine is to be used in conjunction with MatSeqDenseSetPreallocation; 1651 it asserts that the preallocation has a leading dimension (the LDA parameter 1652 of Blas and Lapack fame) larger than M, the first dimension of the matrix. 1653 1654 Level: intermediate 1655 1656 .keywords: dense, matrix, LAPACK, BLAS 1657 1658 .seealso: MatCreate(), MatCreateSeqDense(), MatSeqDenseSetPreallocation() 1659 @*/ 1660 PetscErrorCode MatSeqDenseSetLDA(Mat B,int lda) 1661 { 1662 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1663 PetscFunctionBegin; 1664 if (lda < B->m) SETERRQ2(PETSC_ERR_ARG_SIZ,"LDA %d must be at least matrix dimension %d",lda,B->m); 1665 b->lda = lda; 1666 PetscFunctionReturn(0); 1667 } 1668 1669 /*MC 1670 MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices. 1671 1672 Options Database Keys: 1673 . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions() 1674 1675 Level: beginner 1676 1677 .seealso: MatCreateSeqDense 1678 M*/ 1679 1680 EXTERN_C_BEGIN 1681 #undef __FUNCT__ 1682 #define __FUNCT__ "MatCreate_SeqDense" 1683 PetscErrorCode MatCreate_SeqDense(Mat B) 1684 { 1685 Mat_SeqDense *b; 1686 PetscErrorCode ierr; 1687 int size; 1688 1689 PetscFunctionBegin; 1690 ierr = MPI_Comm_size(B->comm,&size);CHKERRQ(ierr); 1691 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 1692 1693 B->m = B->M = PetscMax(B->m,B->M); 1694 B->n = B->N = PetscMax(B->n,B->N); 1695 1696 ierr = PetscNew(Mat_SeqDense,&b);CHKERRQ(ierr); 1697 ierr = PetscMemzero(b,sizeof(Mat_SeqDense));CHKERRQ(ierr); 1698 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1699 B->factor = 0; 1700 B->mapping = 0; 1701 PetscLogObjectMemory(B,sizeof(struct _p_Mat)); 1702 B->data = (void*)b; 1703 1704 ierr = PetscMapCreateMPI(B->comm,B->m,B->m,&B->rmap);CHKERRQ(ierr); 1705 ierr = PetscMapCreateMPI(B->comm,B->n,B->n,&B->cmap);CHKERRQ(ierr); 1706 1707 b->pivots = 0; 1708 b->roworiented = PETSC_TRUE; 1709 b->v = 0; 1710 b->lda = B->m; 1711 1712 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqDenseSetPreallocation_C", 1713 "MatSeqDenseSetPreallocation_SeqDense", 1714 MatSeqDenseSetPreallocation_SeqDense);CHKERRQ(ierr); 1715 PetscFunctionReturn(0); 1716 } 1717 EXTERN_C_END 1718