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 PetscInt 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 PetscInt 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 PetscInt lda = (PetscInt)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 PetscInt 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(PetscInt)); 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",(PetscInt)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,PetscInt its,PetscInt 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 PetscInt 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 PetscInt _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 PetscInt _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,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **vals) 538 { 539 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 540 PetscScalar *v; 541 PetscErrorCode ierr; 542 PetscInt i; 543 544 PetscFunctionBegin; 545 *ncols = A->n; 546 if (cols) { 547 ierr = PetscMalloc((A->n+1)*sizeof(PetscInt),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,PetscInt row,PetscInt *ncols,PetscInt **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,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],const PetscScalar v[],InsertMode addv) 572 { 573 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 574 PetscInt 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[i],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[j],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[i],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[j],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,PetscInt m,const PetscInt indexm[],PetscInt n,const PetscInt indexn[],PetscScalar v[]) 644 { 645 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 646 PetscInt 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 PetscInt *scols,i,j,nz,header[4]; 671 int fd; 672 PetscMPIInt size; 673 PetscInt *rowlengths = 0,M,N,*cols; 674 PetscScalar *vals,*svals,*v,*w; 675 MPI_Comm comm = ((PetscObject)viewer)->comm; 676 677 PetscFunctionBegin; 678 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 679 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"view must have one processor"); 680 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 681 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 682 if (header[0] != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Not matrix object"); 683 M = header[1]; N = header[2]; nz = header[3]; 684 685 if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense */ 686 ierr = MatCreate(comm,M,N,M,N,A);CHKERRQ(ierr); 687 ierr = MatSetType(*A,type);CHKERRQ(ierr); 688 ierr = MatSeqDenseSetPreallocation(*A,PETSC_NULL);CHKERRQ(ierr); 689 B = *A; 690 a = (Mat_SeqDense*)B->data; 691 v = a->v; 692 /* Allocate some temp space to read in the values and then flip them 693 from row major to column major */ 694 ierr = PetscMalloc((M*N > 0 ? M*N : 1)*sizeof(PetscScalar),&w);CHKERRQ(ierr); 695 /* read in nonzero values */ 696 ierr = PetscBinaryRead(fd,w,M*N,PETSC_SCALAR);CHKERRQ(ierr); 697 /* now flip the values and store them in the matrix*/ 698 for (j=0; j<N; j++) { 699 for (i=0; i<M; i++) { 700 *v++ =w[i*N+j]; 701 } 702 } 703 ierr = PetscFree(w);CHKERRQ(ierr); 704 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 705 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 706 } else { 707 /* read row lengths */ 708 ierr = PetscMalloc((M+1)*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 709 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 710 711 /* create our matrix */ 712 ierr = MatCreate(comm,M,N,M,N,A);CHKERRQ(ierr); 713 ierr = MatSetType(*A,type);CHKERRQ(ierr); 714 ierr = MatSeqDenseSetPreallocation(*A,PETSC_NULL);CHKERRQ(ierr); 715 B = *A; 716 a = (Mat_SeqDense*)B->data; 717 v = a->v; 718 719 /* read column indices and nonzeros */ 720 ierr = PetscMalloc((nz+1)*sizeof(PetscInt),&scols);CHKERRQ(ierr); 721 cols = scols; 722 ierr = PetscBinaryRead(fd,cols,nz,PETSC_INT);CHKERRQ(ierr); 723 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&svals);CHKERRQ(ierr); 724 vals = svals; 725 ierr = PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);CHKERRQ(ierr); 726 727 /* insert into matrix */ 728 for (i=0; i<M; i++) { 729 for (j=0; j<rowlengths[i]; j++) v[i+M*scols[j]] = svals[j]; 730 svals += rowlengths[i]; scols += rowlengths[i]; 731 } 732 ierr = PetscFree(vals);CHKERRQ(ierr); 733 ierr = PetscFree(cols);CHKERRQ(ierr); 734 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 735 736 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 737 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 738 } 739 PetscFunctionReturn(0); 740 } 741 742 #include "petscsys.h" 743 744 #undef __FUNCT__ 745 #define __FUNCT__ "MatView_SeqDense_ASCII" 746 static PetscErrorCode MatView_SeqDense_ASCII(Mat A,PetscViewer viewer) 747 { 748 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 749 PetscErrorCode ierr; 750 PetscInt i,j; 751 char *name; 752 PetscScalar *v; 753 PetscViewerFormat format; 754 755 PetscFunctionBegin; 756 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 757 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 758 if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 759 PetscFunctionReturn(0); /* do nothing for now */ 760 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 761 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 762 for (i=0; i<A->m; i++) { 763 v = a->v + i; 764 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 765 for (j=0; j<A->n; j++) { 766 #if defined(PETSC_USE_COMPLEX) 767 if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) { 768 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i) ",j,PetscRealPart(*v),PetscImaginaryPart(*v));CHKERRQ(ierr); 769 } else if (PetscRealPart(*v)) { 770 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,PetscRealPart(*v));CHKERRQ(ierr); 771 } 772 #else 773 if (*v) { 774 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",j,*v);CHKERRQ(ierr); 775 } 776 #endif 777 v += a->lda; 778 } 779 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 780 } 781 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 782 } else { 783 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 784 #if defined(PETSC_USE_COMPLEX) 785 PetscTruth allreal = PETSC_TRUE; 786 /* determine if matrix has all real values */ 787 v = a->v; 788 for (i=0; i<A->m*A->n; i++) { 789 if (PetscImaginaryPart(v[i])) { allreal = PETSC_FALSE; break ;} 790 } 791 #endif 792 if (format == PETSC_VIEWER_ASCII_MATLAB) { 793 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 794 ierr = PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",A->m,A->n);CHKERRQ(ierr); 795 ierr = PetscViewerASCIIPrintf(viewer,"%s = zeros(%D,%D);\n",name,A->m,A->n);CHKERRQ(ierr); 796 ierr = PetscViewerASCIIPrintf(viewer,"%s = [\n",name);CHKERRQ(ierr); 797 } 798 799 for (i=0; i<A->m; i++) { 800 v = a->v + i; 801 for (j=0; j<A->n; j++) { 802 #if defined(PETSC_USE_COMPLEX) 803 if (allreal) { 804 ierr = PetscViewerASCIIPrintf(viewer,"%6.4e ",PetscRealPart(*v));CHKERRQ(ierr); 805 } else { 806 ierr = PetscViewerASCIIPrintf(viewer,"%6.4e + %6.4e i ",PetscRealPart(*v),PetscImaginaryPart(*v));CHKERRQ(ierr); 807 } 808 #else 809 ierr = PetscViewerASCIIPrintf(viewer,"%6.4e ",*v);CHKERRQ(ierr); 810 #endif 811 v += a->lda; 812 } 813 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 814 } 815 if (format == PETSC_VIEWER_ASCII_MATLAB) { 816 ierr = PetscViewerASCIIPrintf(viewer,"];\n");CHKERRQ(ierr); 817 } 818 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 819 } 820 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 821 PetscFunctionReturn(0); 822 } 823 824 #undef __FUNCT__ 825 #define __FUNCT__ "MatView_SeqDense_Binary" 826 static PetscErrorCode MatView_SeqDense_Binary(Mat A,PetscViewer viewer) 827 { 828 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 829 PetscErrorCode ierr; 830 int fd; 831 PetscInt ict,j,n = A->n,m = A->m,i,*col_lens,nz = m*n; 832 PetscScalar *v,*anonz,*vals; 833 PetscViewerFormat format; 834 835 PetscFunctionBegin; 836 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 837 838 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 839 if (format == PETSC_VIEWER_BINARY_NATIVE) { 840 /* store the matrix as a dense matrix */ 841 ierr = PetscMalloc(4*sizeof(PetscInt),&col_lens);CHKERRQ(ierr); 842 col_lens[0] = MAT_FILE_COOKIE; 843 col_lens[1] = m; 844 col_lens[2] = n; 845 col_lens[3] = MATRIX_BINARY_FORMAT_DENSE; 846 ierr = PetscBinaryWrite(fd,col_lens,4,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 847 ierr = PetscFree(col_lens);CHKERRQ(ierr); 848 849 /* write out matrix, by rows */ 850 ierr = PetscMalloc((m*n+1)*sizeof(PetscScalar),&vals);CHKERRQ(ierr); 851 v = a->v; 852 for (i=0; i<m; i++) { 853 for (j=0; j<n; j++) { 854 vals[i + j*m] = *v++; 855 } 856 } 857 ierr = PetscBinaryWrite(fd,vals,n*m,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 858 ierr = PetscFree(vals);CHKERRQ(ierr); 859 } else { 860 ierr = PetscMalloc((4+nz)*sizeof(PetscInt),&col_lens);CHKERRQ(ierr); 861 col_lens[0] = MAT_FILE_COOKIE; 862 col_lens[1] = m; 863 col_lens[2] = n; 864 col_lens[3] = nz; 865 866 /* store lengths of each row and write (including header) to file */ 867 for (i=0; i<m; i++) col_lens[4+i] = n; 868 ierr = PetscBinaryWrite(fd,col_lens,4+m,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 869 870 /* Possibly should write in smaller increments, not whole matrix at once? */ 871 /* store column indices (zero start index) */ 872 ict = 0; 873 for (i=0; i<m; i++) { 874 for (j=0; j<n; j++) col_lens[ict++] = j; 875 } 876 ierr = PetscBinaryWrite(fd,col_lens,nz,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 877 ierr = PetscFree(col_lens);CHKERRQ(ierr); 878 879 /* store nonzero values */ 880 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&anonz);CHKERRQ(ierr); 881 ict = 0; 882 for (i=0; i<m; i++) { 883 v = a->v + i; 884 for (j=0; j<n; j++) { 885 anonz[ict++] = *v; v += a->lda; 886 } 887 } 888 ierr = PetscBinaryWrite(fd,anonz,nz,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 889 ierr = PetscFree(anonz);CHKERRQ(ierr); 890 } 891 PetscFunctionReturn(0); 892 } 893 894 #undef __FUNCT__ 895 #define __FUNCT__ "MatView_SeqDense_Draw_Zoom" 896 PetscErrorCode MatView_SeqDense_Draw_Zoom(PetscDraw draw,void *Aa) 897 { 898 Mat A = (Mat) Aa; 899 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 900 PetscErrorCode ierr; 901 PetscInt m = A->m,n = A->n,color,i,j; 902 PetscScalar *v = a->v; 903 PetscViewer viewer; 904 PetscDraw popup; 905 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r,scale,maxv = 0.0; 906 PetscViewerFormat format; 907 908 PetscFunctionBegin; 909 910 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 911 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 912 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 913 914 /* Loop over matrix elements drawing boxes */ 915 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 916 /* Blue for negative and Red for positive */ 917 color = PETSC_DRAW_BLUE; 918 for(j = 0; j < n; j++) { 919 x_l = j; 920 x_r = x_l + 1.0; 921 for(i = 0; i < m; i++) { 922 y_l = m - i - 1.0; 923 y_r = y_l + 1.0; 924 #if defined(PETSC_USE_COMPLEX) 925 if (PetscRealPart(v[j*m+i]) > 0.) { 926 color = PETSC_DRAW_RED; 927 } else if (PetscRealPart(v[j*m+i]) < 0.) { 928 color = PETSC_DRAW_BLUE; 929 } else { 930 continue; 931 } 932 #else 933 if (v[j*m+i] > 0.) { 934 color = PETSC_DRAW_RED; 935 } else if (v[j*m+i] < 0.) { 936 color = PETSC_DRAW_BLUE; 937 } else { 938 continue; 939 } 940 #endif 941 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 942 } 943 } 944 } else { 945 /* use contour shading to indicate magnitude of values */ 946 /* first determine max of all nonzero values */ 947 for(i = 0; i < m*n; i++) { 948 if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]); 949 } 950 scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/maxv; 951 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 952 if (popup) {ierr = PetscDrawScalePopup(popup,0.0,maxv);CHKERRQ(ierr);} 953 for(j = 0; j < n; j++) { 954 x_l = j; 955 x_r = x_l + 1.0; 956 for(i = 0; i < m; i++) { 957 y_l = m - i - 1.0; 958 y_r = y_l + 1.0; 959 color = PETSC_DRAW_BASIC_COLORS + (int)(scale*PetscAbsScalar(v[j*m+i])); 960 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 961 } 962 } 963 } 964 PetscFunctionReturn(0); 965 } 966 967 #undef __FUNCT__ 968 #define __FUNCT__ "MatView_SeqDense_Draw" 969 PetscErrorCode MatView_SeqDense_Draw(Mat A,PetscViewer viewer) 970 { 971 PetscDraw draw; 972 PetscTruth isnull; 973 PetscReal xr,yr,xl,yl,h,w; 974 PetscErrorCode ierr; 975 976 PetscFunctionBegin; 977 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 978 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 979 if (isnull == PETSC_TRUE) PetscFunctionReturn(0); 980 981 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 982 xr = A->n; yr = A->m; h = yr/10.0; w = xr/10.0; 983 xr += w; yr += h; xl = -w; yl = -h; 984 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 985 ierr = PetscDrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);CHKERRQ(ierr); 986 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);CHKERRQ(ierr); 987 PetscFunctionReturn(0); 988 } 989 990 #undef __FUNCT__ 991 #define __FUNCT__ "MatView_SeqDense" 992 PetscErrorCode MatView_SeqDense(Mat A,PetscViewer viewer) 993 { 994 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 995 PetscErrorCode ierr; 996 PetscTruth issocket,iascii,isbinary,isdraw; 997 998 PetscFunctionBegin; 999 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_SOCKET,&issocket);CHKERRQ(ierr); 1000 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 1001 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_BINARY,&isbinary);CHKERRQ(ierr); 1002 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_DRAW,&isdraw);CHKERRQ(ierr); 1003 1004 if (issocket) { 1005 if (a->lda>A->m) SETERRQ(PETSC_ERR_SUP,"Case can not handle LDA"); 1006 ierr = PetscViewerSocketPutScalar(viewer,A->m,A->n,a->v);CHKERRQ(ierr); 1007 } else if (iascii) { 1008 ierr = MatView_SeqDense_ASCII(A,viewer);CHKERRQ(ierr); 1009 } else if (isbinary) { 1010 ierr = MatView_SeqDense_Binary(A,viewer);CHKERRQ(ierr); 1011 } else if (isdraw) { 1012 ierr = MatView_SeqDense_Draw(A,viewer);CHKERRQ(ierr); 1013 } else { 1014 SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported by dense matrix",((PetscObject)viewer)->type_name); 1015 } 1016 PetscFunctionReturn(0); 1017 } 1018 1019 #undef __FUNCT__ 1020 #define __FUNCT__ "MatDestroy_SeqDense" 1021 PetscErrorCode MatDestroy_SeqDense(Mat mat) 1022 { 1023 Mat_SeqDense *l = (Mat_SeqDense*)mat->data; 1024 PetscErrorCode ierr; 1025 1026 PetscFunctionBegin; 1027 #if defined(PETSC_USE_LOG) 1028 PetscLogObjectState((PetscObject)mat,"Rows %D Cols %D",mat->m,mat->n); 1029 #endif 1030 if (l->pivots) {ierr = PetscFree(l->pivots);CHKERRQ(ierr);} 1031 if (!l->user_alloc) {ierr = PetscFree(l->v);CHKERRQ(ierr);} 1032 ierr = PetscFree(l);CHKERRQ(ierr); 1033 ierr = PetscObjectComposeFunctionDynamic((PetscObject)mat,"MatSeqDenseSetPreallocation_C","",PETSC_NULL);CHKERRQ(ierr); 1034 PetscFunctionReturn(0); 1035 } 1036 1037 #undef __FUNCT__ 1038 #define __FUNCT__ "MatTranspose_SeqDense" 1039 PetscErrorCode MatTranspose_SeqDense(Mat A,Mat *matout) 1040 { 1041 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1042 PetscErrorCode ierr; 1043 PetscInt k,j,m,n,M; 1044 PetscScalar *v,tmp; 1045 1046 PetscFunctionBegin; 1047 v = mat->v; m = A->m; M = mat->lda; n = A->n; 1048 if (!matout) { /* in place transpose */ 1049 if (m != n) { 1050 SETERRQ(PETSC_ERR_SUP,"Can not transpose non-square matrix in place"); 1051 } else { 1052 for (j=0; j<m; j++) { 1053 for (k=0; k<j; k++) { 1054 tmp = v[j + k*M]; 1055 v[j + k*M] = v[k + j*M]; 1056 v[k + j*M] = tmp; 1057 } 1058 } 1059 } 1060 } else { /* out-of-place transpose */ 1061 Mat tmat; 1062 Mat_SeqDense *tmatd; 1063 PetscScalar *v2; 1064 1065 ierr = MatCreate(A->comm,A->n,A->m,A->n,A->m,&tmat);CHKERRQ(ierr); 1066 ierr = MatSetType(tmat,A->type_name);CHKERRQ(ierr); 1067 ierr = MatSeqDenseSetPreallocation(tmat,PETSC_NULL);CHKERRQ(ierr); 1068 tmatd = (Mat_SeqDense*)tmat->data; 1069 v = mat->v; v2 = tmatd->v; 1070 for (j=0; j<n; j++) { 1071 for (k=0; k<m; k++) v2[j + k*n] = v[k + j*M]; 1072 } 1073 ierr = MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1074 ierr = MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1075 *matout = tmat; 1076 } 1077 PetscFunctionReturn(0); 1078 } 1079 1080 #undef __FUNCT__ 1081 #define __FUNCT__ "MatEqual_SeqDense" 1082 PetscErrorCode MatEqual_SeqDense(Mat A1,Mat A2,PetscTruth *flg) 1083 { 1084 Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data; 1085 Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data; 1086 PetscInt i,j; 1087 PetscScalar *v1 = mat1->v,*v2 = mat2->v; 1088 1089 PetscFunctionBegin; 1090 if (A1->m != A2->m) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1091 if (A1->n != A2->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1092 for (i=0; i<A1->m; i++) { 1093 v1 = mat1->v+i; v2 = mat2->v+i; 1094 for (j=0; j<A1->n; j++) { 1095 if (*v1 != *v2) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 1096 v1 += mat1->lda; v2 += mat2->lda; 1097 } 1098 } 1099 *flg = PETSC_TRUE; 1100 PetscFunctionReturn(0); 1101 } 1102 1103 #undef __FUNCT__ 1104 #define __FUNCT__ "MatGetDiagonal_SeqDense" 1105 PetscErrorCode MatGetDiagonal_SeqDense(Mat A,Vec v) 1106 { 1107 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1108 PetscErrorCode ierr; 1109 PetscInt i,n,len; 1110 PetscScalar *x,zero = 0.0; 1111 1112 PetscFunctionBegin; 1113 ierr = VecSet(&zero,v);CHKERRQ(ierr); 1114 ierr = VecGetSize(v,&n);CHKERRQ(ierr); 1115 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1116 len = PetscMin(A->m,A->n); 1117 if (n != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming mat and vec"); 1118 for (i=0; i<len; i++) { 1119 x[i] = mat->v[i*mat->lda + i]; 1120 } 1121 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1122 PetscFunctionReturn(0); 1123 } 1124 1125 #undef __FUNCT__ 1126 #define __FUNCT__ "MatDiagonalScale_SeqDense" 1127 PetscErrorCode MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr) 1128 { 1129 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1130 PetscScalar *l,*r,x,*v; 1131 PetscErrorCode ierr; 1132 PetscInt i,j,m = A->m,n = A->n; 1133 1134 PetscFunctionBegin; 1135 if (ll) { 1136 ierr = VecGetSize(ll,&m);CHKERRQ(ierr); 1137 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1138 if (m != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Left scaling vec wrong size"); 1139 for (i=0; i<m; i++) { 1140 x = l[i]; 1141 v = mat->v + i; 1142 for (j=0; j<n; j++) { (*v) *= x; v+= m;} 1143 } 1144 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1145 PetscLogFlops(n*m); 1146 } 1147 if (rr) { 1148 ierr = VecGetSize(rr,&n);CHKERRQ(ierr); 1149 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1150 if (n != A->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Right scaling vec wrong size"); 1151 for (i=0; i<n; i++) { 1152 x = r[i]; 1153 v = mat->v + i*m; 1154 for (j=0; j<m; j++) { (*v++) *= x;} 1155 } 1156 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1157 PetscLogFlops(n*m); 1158 } 1159 PetscFunctionReturn(0); 1160 } 1161 1162 #undef __FUNCT__ 1163 #define __FUNCT__ "MatNorm_SeqDense" 1164 PetscErrorCode MatNorm_SeqDense(Mat A,NormType type,PetscReal *nrm) 1165 { 1166 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1167 PetscScalar *v = mat->v; 1168 PetscReal sum = 0.0; 1169 PetscInt lda=mat->lda,m=A->m,i,j; 1170 1171 PetscFunctionBegin; 1172 if (type == NORM_FROBENIUS) { 1173 if (lda>m) { 1174 for (j=0; j<A->n; j++) { 1175 v = mat->v+j*lda; 1176 for (i=0; i<m; i++) { 1177 #if defined(PETSC_USE_COMPLEX) 1178 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1179 #else 1180 sum += (*v)*(*v); v++; 1181 #endif 1182 } 1183 } 1184 } else { 1185 for (i=0; i<A->n*A->m; i++) { 1186 #if defined(PETSC_USE_COMPLEX) 1187 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1188 #else 1189 sum += (*v)*(*v); v++; 1190 #endif 1191 } 1192 } 1193 *nrm = sqrt(sum); 1194 PetscLogFlops(2*A->n*A->m); 1195 } else if (type == NORM_1) { 1196 *nrm = 0.0; 1197 for (j=0; j<A->n; j++) { 1198 v = mat->v + j*mat->lda; 1199 sum = 0.0; 1200 for (i=0; i<A->m; i++) { 1201 sum += PetscAbsScalar(*v); v++; 1202 } 1203 if (sum > *nrm) *nrm = sum; 1204 } 1205 PetscLogFlops(A->n*A->m); 1206 } else if (type == NORM_INFINITY) { 1207 *nrm = 0.0; 1208 for (j=0; j<A->m; j++) { 1209 v = mat->v + j; 1210 sum = 0.0; 1211 for (i=0; i<A->n; i++) { 1212 sum += PetscAbsScalar(*v); v += mat->lda; 1213 } 1214 if (sum > *nrm) *nrm = sum; 1215 } 1216 PetscLogFlops(A->n*A->m); 1217 } else { 1218 SETERRQ(PETSC_ERR_SUP,"No two norm"); 1219 } 1220 PetscFunctionReturn(0); 1221 } 1222 1223 #undef __FUNCT__ 1224 #define __FUNCT__ "MatSetOption_SeqDense" 1225 PetscErrorCode MatSetOption_SeqDense(Mat A,MatOption op) 1226 { 1227 Mat_SeqDense *aij = (Mat_SeqDense*)A->data; 1228 1229 PetscFunctionBegin; 1230 switch (op) { 1231 case MAT_ROW_ORIENTED: 1232 aij->roworiented = PETSC_TRUE; 1233 break; 1234 case MAT_COLUMN_ORIENTED: 1235 aij->roworiented = PETSC_FALSE; 1236 break; 1237 case MAT_ROWS_SORTED: 1238 case MAT_ROWS_UNSORTED: 1239 case MAT_COLUMNS_SORTED: 1240 case MAT_COLUMNS_UNSORTED: 1241 case MAT_NO_NEW_NONZERO_LOCATIONS: 1242 case MAT_YES_NEW_NONZERO_LOCATIONS: 1243 case MAT_NEW_NONZERO_LOCATION_ERR: 1244 case MAT_NO_NEW_DIAGONALS: 1245 case MAT_YES_NEW_DIAGONALS: 1246 case MAT_IGNORE_OFF_PROC_ENTRIES: 1247 case MAT_USE_HASH_TABLE: 1248 PetscLogInfo(A,"MatSetOption_SeqDense:Option ignored\n"); 1249 break; 1250 case MAT_SYMMETRIC: 1251 case MAT_STRUCTURALLY_SYMMETRIC: 1252 case MAT_NOT_SYMMETRIC: 1253 case MAT_NOT_STRUCTURALLY_SYMMETRIC: 1254 case MAT_HERMITIAN: 1255 case MAT_NOT_HERMITIAN: 1256 case MAT_SYMMETRY_ETERNAL: 1257 case MAT_NOT_SYMMETRY_ETERNAL: 1258 break; 1259 default: 1260 SETERRQ(PETSC_ERR_SUP,"unknown option"); 1261 } 1262 PetscFunctionReturn(0); 1263 } 1264 1265 #undef __FUNCT__ 1266 #define __FUNCT__ "MatZeroEntries_SeqDense" 1267 PetscErrorCode MatZeroEntries_SeqDense(Mat A) 1268 { 1269 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1270 PetscErrorCode ierr; 1271 PetscInt lda=l->lda,m=A->m,j; 1272 1273 PetscFunctionBegin; 1274 if (lda>m) { 1275 for (j=0; j<A->n; j++) { 1276 ierr = PetscMemzero(l->v+j*lda,m*sizeof(PetscScalar));CHKERRQ(ierr); 1277 } 1278 } else { 1279 ierr = PetscMemzero(l->v,A->m*A->n*sizeof(PetscScalar));CHKERRQ(ierr); 1280 } 1281 PetscFunctionReturn(0); 1282 } 1283 1284 #undef __FUNCT__ 1285 #define __FUNCT__ "MatZeroRows_SeqDense" 1286 PetscErrorCode MatZeroRows_SeqDense(Mat A,IS is,const PetscScalar *diag) 1287 { 1288 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1289 PetscErrorCode ierr; 1290 PetscInt n = A->n,i,j,N,*rows; 1291 PetscScalar *slot; 1292 1293 PetscFunctionBegin; 1294 ierr = ISGetLocalSize(is,&N);CHKERRQ(ierr); 1295 ierr = ISGetIndices(is,&rows);CHKERRQ(ierr); 1296 for (i=0; i<N; i++) { 1297 slot = l->v + rows[i]; 1298 for (j=0; j<n; j++) { *slot = 0.0; slot += n;} 1299 } 1300 if (diag) { 1301 for (i=0; i<N; i++) { 1302 slot = l->v + (n+1)*rows[i]; 1303 *slot = *diag; 1304 } 1305 } 1306 ISRestoreIndices(is,&rows); 1307 PetscFunctionReturn(0); 1308 } 1309 1310 #undef __FUNCT__ 1311 #define __FUNCT__ "MatGetArray_SeqDense" 1312 PetscErrorCode MatGetArray_SeqDense(Mat A,PetscScalar *array[]) 1313 { 1314 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1315 1316 PetscFunctionBegin; 1317 *array = mat->v; 1318 PetscFunctionReturn(0); 1319 } 1320 1321 #undef __FUNCT__ 1322 #define __FUNCT__ "MatRestoreArray_SeqDense" 1323 PetscErrorCode MatRestoreArray_SeqDense(Mat A,PetscScalar *array[]) 1324 { 1325 PetscFunctionBegin; 1326 *array = 0; /* user cannot accidently use the array later */ 1327 PetscFunctionReturn(0); 1328 } 1329 1330 #undef __FUNCT__ 1331 #define __FUNCT__ "MatGetSubMatrix_SeqDense" 1332 static PetscErrorCode MatGetSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,PetscInt cs,MatReuse scall,Mat *B) 1333 { 1334 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1335 PetscErrorCode ierr; 1336 PetscInt i,j,m = A->m,*irow,*icol,nrows,ncols; 1337 PetscScalar *av,*bv,*v = mat->v; 1338 Mat newmat; 1339 1340 PetscFunctionBegin; 1341 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1342 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1343 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1344 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1345 1346 /* Check submatrixcall */ 1347 if (scall == MAT_REUSE_MATRIX) { 1348 PetscInt n_cols,n_rows; 1349 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1350 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 1351 newmat = *B; 1352 } else { 1353 /* Create and fill new matrix */ 1354 ierr = MatCreate(A->comm,nrows,ncols,nrows,ncols,&newmat);CHKERRQ(ierr); 1355 ierr = MatSetType(newmat,A->type_name);CHKERRQ(ierr); 1356 ierr = MatSeqDenseSetPreallocation(newmat,PETSC_NULL);CHKERRQ(ierr); 1357 } 1358 1359 /* Now extract the data pointers and do the copy,column at a time */ 1360 bv = ((Mat_SeqDense*)newmat->data)->v; 1361 1362 for (i=0; i<ncols; i++) { 1363 av = v + m*icol[i]; 1364 for (j=0; j<nrows; j++) { 1365 *bv++ = av[irow[j]]; 1366 } 1367 } 1368 1369 /* Assemble the matrices so that the correct flags are set */ 1370 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1371 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1372 1373 /* Free work space */ 1374 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1375 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1376 *B = newmat; 1377 PetscFunctionReturn(0); 1378 } 1379 1380 #undef __FUNCT__ 1381 #define __FUNCT__ "MatGetSubMatrices_SeqDense" 1382 PetscErrorCode MatGetSubMatrices_SeqDense(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1383 { 1384 PetscErrorCode ierr; 1385 PetscInt i; 1386 1387 PetscFunctionBegin; 1388 if (scall == MAT_INITIAL_MATRIX) { 1389 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1390 } 1391 1392 for (i=0; i<n; i++) { 1393 ierr = MatGetSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1394 } 1395 PetscFunctionReturn(0); 1396 } 1397 1398 #undef __FUNCT__ 1399 #define __FUNCT__ "MatCopy_SeqDense" 1400 PetscErrorCode MatCopy_SeqDense(Mat A,Mat B,MatStructure str) 1401 { 1402 Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense *)B->data; 1403 PetscErrorCode ierr; 1404 PetscInt lda1=a->lda,lda2=b->lda, m=A->m,n=A->n, j; 1405 1406 PetscFunctionBegin; 1407 /* If the two matrices don't have the same copy implementation, they aren't compatible for fast copy. */ 1408 if (A->ops->copy != B->ops->copy) { 1409 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1410 PetscFunctionReturn(0); 1411 } 1412 if (m != B->m || n != B->n) SETERRQ(PETSC_ERR_ARG_SIZ,"size(B) != size(A)"); 1413 if (lda1>m || lda2>m) { 1414 for (j=0; j<n; j++) { 1415 ierr = PetscMemcpy(b->v+j*lda2,a->v+j*lda1,m*sizeof(PetscScalar));CHKERRQ(ierr); 1416 } 1417 } else { 1418 ierr = PetscMemcpy(b->v,a->v,A->m*A->n*sizeof(PetscScalar));CHKERRQ(ierr); 1419 } 1420 PetscFunctionReturn(0); 1421 } 1422 1423 #undef __FUNCT__ 1424 #define __FUNCT__ "MatSetUpPreallocation_SeqDense" 1425 PetscErrorCode MatSetUpPreallocation_SeqDense(Mat A) 1426 { 1427 PetscErrorCode ierr; 1428 1429 PetscFunctionBegin; 1430 ierr = MatSeqDenseSetPreallocation(A,0);CHKERRQ(ierr); 1431 PetscFunctionReturn(0); 1432 } 1433 1434 /* -------------------------------------------------------------------*/ 1435 static struct _MatOps MatOps_Values = {MatSetValues_SeqDense, 1436 MatGetRow_SeqDense, 1437 MatRestoreRow_SeqDense, 1438 MatMult_SeqDense, 1439 /* 4*/ MatMultAdd_SeqDense, 1440 MatMultTranspose_SeqDense, 1441 MatMultTransposeAdd_SeqDense, 1442 MatSolve_SeqDense, 1443 MatSolveAdd_SeqDense, 1444 MatSolveTranspose_SeqDense, 1445 /*10*/ MatSolveTransposeAdd_SeqDense, 1446 MatLUFactor_SeqDense, 1447 MatCholeskyFactor_SeqDense, 1448 MatRelax_SeqDense, 1449 MatTranspose_SeqDense, 1450 /*15*/ MatGetInfo_SeqDense, 1451 MatEqual_SeqDense, 1452 MatGetDiagonal_SeqDense, 1453 MatDiagonalScale_SeqDense, 1454 MatNorm_SeqDense, 1455 /*20*/ 0, 1456 0, 1457 0, 1458 MatSetOption_SeqDense, 1459 MatZeroEntries_SeqDense, 1460 /*25*/ MatZeroRows_SeqDense, 1461 MatLUFactorSymbolic_SeqDense, 1462 MatLUFactorNumeric_SeqDense, 1463 MatCholeskyFactorSymbolic_SeqDense, 1464 MatCholeskyFactorNumeric_SeqDense, 1465 /*30*/ MatSetUpPreallocation_SeqDense, 1466 0, 1467 0, 1468 MatGetArray_SeqDense, 1469 MatRestoreArray_SeqDense, 1470 /*35*/ MatDuplicate_SeqDense, 1471 0, 1472 0, 1473 0, 1474 0, 1475 /*40*/ MatAXPY_SeqDense, 1476 MatGetSubMatrices_SeqDense, 1477 0, 1478 MatGetValues_SeqDense, 1479 MatCopy_SeqDense, 1480 /*45*/ 0, 1481 MatScale_SeqDense, 1482 0, 1483 0, 1484 0, 1485 /*50*/ 0, 1486 0, 1487 0, 1488 0, 1489 0, 1490 /*55*/ 0, 1491 0, 1492 0, 1493 0, 1494 0, 1495 /*60*/ 0, 1496 MatDestroy_SeqDense, 1497 MatView_SeqDense, 1498 MatGetPetscMaps_Petsc, 1499 0, 1500 /*65*/ 0, 1501 0, 1502 0, 1503 0, 1504 0, 1505 /*70*/ 0, 1506 0, 1507 0, 1508 0, 1509 0, 1510 /*75*/ 0, 1511 0, 1512 0, 1513 0, 1514 0, 1515 /*80*/ 0, 1516 0, 1517 0, 1518 0, 1519 /*84*/ MatLoad_SeqDense, 1520 0, 1521 0, 1522 0, 1523 0, 1524 0, 1525 /*90*/ 0, 1526 0, 1527 0, 1528 0, 1529 0, 1530 /*95*/ 0, 1531 0, 1532 0, 1533 0}; 1534 1535 #undef __FUNCT__ 1536 #define __FUNCT__ "MatCreateSeqDense" 1537 /*@C 1538 MatCreateSeqDense - Creates a sequential dense matrix that 1539 is stored in column major order (the usual Fortran 77 manner). Many 1540 of the matrix operations use the BLAS and LAPACK routines. 1541 1542 Collective on MPI_Comm 1543 1544 Input Parameters: 1545 + comm - MPI communicator, set to PETSC_COMM_SELF 1546 . m - number of rows 1547 . n - number of columns 1548 - data - optional location of matrix data. Set data=PETSC_NULL for PETSc 1549 to control all matrix memory allocation. 1550 1551 Output Parameter: 1552 . A - the matrix 1553 1554 Notes: 1555 The data input variable is intended primarily for Fortran programmers 1556 who wish to allocate their own matrix memory space. Most users should 1557 set data=PETSC_NULL. 1558 1559 Level: intermediate 1560 1561 .keywords: dense, matrix, LAPACK, BLAS 1562 1563 .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues() 1564 @*/ 1565 PetscErrorCode MatCreateSeqDense(MPI_Comm comm,PetscInt m,PetscInt n,PetscScalar *data,Mat *A) 1566 { 1567 PetscErrorCode ierr; 1568 1569 PetscFunctionBegin; 1570 ierr = MatCreate(comm,m,n,m,n,A);CHKERRQ(ierr); 1571 ierr = MatSetType(*A,MATSEQDENSE);CHKERRQ(ierr); 1572 ierr = MatSeqDenseSetPreallocation(*A,data);CHKERRQ(ierr); 1573 PetscFunctionReturn(0); 1574 } 1575 1576 #undef __FUNCT__ 1577 #define __FUNCT__ "MatSeqDensePreallocation" 1578 /*@C 1579 MatSeqDenseSetPreallocation - Sets the array used for storing the matrix elements 1580 1581 Collective on MPI_Comm 1582 1583 Input Parameters: 1584 + A - the matrix 1585 - data - the array (or PETSC_NULL) 1586 1587 Notes: 1588 The data input variable is intended primarily for Fortran programmers 1589 who wish to allocate their own matrix memory space. Most users should 1590 set data=PETSC_NULL. 1591 1592 Level: intermediate 1593 1594 .keywords: dense, matrix, LAPACK, BLAS 1595 1596 .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues() 1597 @*/ 1598 PetscErrorCode MatSeqDenseSetPreallocation(Mat B,PetscScalar data[]) 1599 { 1600 PetscErrorCode ierr,(*f)(Mat,PetscScalar[]); 1601 1602 PetscFunctionBegin; 1603 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqDenseSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 1604 if (f) { 1605 ierr = (*f)(B,data);CHKERRQ(ierr); 1606 } 1607 PetscFunctionReturn(0); 1608 } 1609 1610 EXTERN_C_BEGIN 1611 #undef __FUNCT__ 1612 #define __FUNCT__ "MatSeqDensePreallocation_SeqDense" 1613 PetscErrorCode MatSeqDenseSetPreallocation_SeqDense(Mat B,PetscScalar *data) 1614 { 1615 Mat_SeqDense *b; 1616 PetscErrorCode ierr; 1617 1618 PetscFunctionBegin; 1619 B->preallocated = PETSC_TRUE; 1620 b = (Mat_SeqDense*)B->data; 1621 if (!data) { 1622 ierr = PetscMalloc((B->m*B->n+1)*sizeof(PetscScalar),&b->v);CHKERRQ(ierr); 1623 ierr = PetscMemzero(b->v,B->m*B->n*sizeof(PetscScalar));CHKERRQ(ierr); 1624 b->user_alloc = PETSC_FALSE; 1625 PetscLogObjectMemory(B,B->n*B->m*sizeof(PetscScalar)); 1626 } else { /* user-allocated storage */ 1627 b->v = data; 1628 b->user_alloc = PETSC_TRUE; 1629 } 1630 PetscFunctionReturn(0); 1631 } 1632 EXTERN_C_END 1633 1634 #undef __FUNCT__ 1635 #define __FUNCT__ "MatSeqDenseSetLDA" 1636 /*@C 1637 MatSeqDenseSetLDA - Declare the leading dimension of the user-provided array 1638 1639 Input parameter: 1640 + A - the matrix 1641 - lda - the leading dimension 1642 1643 Notes: 1644 This routine is to be used in conjunction with MatSeqDenseSetPreallocation; 1645 it asserts that the preallocation has a leading dimension (the LDA parameter 1646 of Blas and Lapack fame) larger than M, the first dimension of the matrix. 1647 1648 Level: intermediate 1649 1650 .keywords: dense, matrix, LAPACK, BLAS 1651 1652 .seealso: MatCreate(), MatCreateSeqDense(), MatSeqDenseSetPreallocation() 1653 @*/ 1654 PetscErrorCode MatSeqDenseSetLDA(Mat B,PetscInt lda) 1655 { 1656 Mat_SeqDense *b = (Mat_SeqDense*)B->data; 1657 PetscFunctionBegin; 1658 if (lda < B->m) SETERRQ2(PETSC_ERR_ARG_SIZ,"LDA %D must be at least matrix dimension %D",lda,B->m); 1659 b->lda = lda; 1660 PetscFunctionReturn(0); 1661 } 1662 1663 /*MC 1664 MATSEQDENSE - MATSEQDENSE = "seqdense" - A matrix type to be used for sequential dense matrices. 1665 1666 Options Database Keys: 1667 . -mat_type seqdense - sets the matrix type to "seqdense" during a call to MatSetFromOptions() 1668 1669 Level: beginner 1670 1671 .seealso: MatCreateSeqDense 1672 M*/ 1673 1674 EXTERN_C_BEGIN 1675 #undef __FUNCT__ 1676 #define __FUNCT__ "MatCreate_SeqDense" 1677 PetscErrorCode MatCreate_SeqDense(Mat B) 1678 { 1679 Mat_SeqDense *b; 1680 PetscErrorCode ierr; 1681 PetscMPIInt size; 1682 1683 PetscFunctionBegin; 1684 ierr = MPI_Comm_size(B->comm,&size);CHKERRQ(ierr); 1685 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 1686 1687 B->m = B->M = PetscMax(B->m,B->M); 1688 B->n = B->N = PetscMax(B->n,B->N); 1689 1690 ierr = PetscNew(Mat_SeqDense,&b);CHKERRQ(ierr); 1691 ierr = PetscMemzero(b,sizeof(Mat_SeqDense));CHKERRQ(ierr); 1692 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1693 B->factor = 0; 1694 B->mapping = 0; 1695 PetscLogObjectMemory(B,sizeof(struct _p_Mat)); 1696 B->data = (void*)b; 1697 1698 ierr = PetscMapCreateMPI(B->comm,B->m,B->m,&B->rmap);CHKERRQ(ierr); 1699 ierr = PetscMapCreateMPI(B->comm,B->n,B->n,&B->cmap);CHKERRQ(ierr); 1700 1701 b->pivots = 0; 1702 b->roworiented = PETSC_TRUE; 1703 b->v = 0; 1704 b->lda = B->m; 1705 1706 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqDenseSetPreallocation_C", 1707 "MatSeqDenseSetPreallocation_SeqDense", 1708 MatSeqDenseSetPreallocation_SeqDense);CHKERRQ(ierr); 1709 PetscFunctionReturn(0); 1710 } 1711 EXTERN_C_END 1712