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