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