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