1 /*$Id: dense.c,v 1.184 2000/04/09 04:35: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 "pinclude/blaslapack.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 inA) 54 { 55 Mat_SeqDense *a = (Mat_SeqDense*)inA->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,PetscReal f) 71 { 72 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 73 int info; 74 75 PetscFunctionBegin; 76 if (!mat->pivots) { 77 mat->pivots = (int*)PetscMalloc((mat->m+1)*sizeof(int));CHKPTRQ(mat->pivots); 78 PLogObjectMemory(A,mat->m*sizeof(int)); 79 } 80 A->factor = FACTOR_LU; 81 if (!mat->m || !mat->n) PetscFunctionReturn(0); 82 LAgetrf_(&mat->m,&mat->n,mat->v,&mat->m,mat->pivots,&info); 83 if (info<0) SETERRQ(PETSC_ERR_LIB,0,"Bad argument to LU factorization"); 84 if (info>0) SETERRQ(PETSC_ERR_MAT_LU_ZRPVT,0,"Bad LU factorization"); 85 PLogFlops((2*mat->n*mat->n*mat->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,mat->m,mat->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,mat->m*mat->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,PetscReal f,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,mat->m*mat->n*sizeof(Scalar));CHKERRQ(ierr); 129 (*fact)->factor = 0; 130 ierr = MatLUFactor(*fact,0,0,1.0);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,-mat->m*sizeof(int)); 167 mat->pivots = 0; 168 } 169 if (!mat->m || !mat->n) PetscFunctionReturn(0); 170 LApotrf_("L",&mat->n,mat->v,&mat->m,&info); 171 if (info) SETERRQ1(PETSC_ERR_MAT_CH_ZRPVT,0,"Bad factorization: zero pivot in row %d",info-1); 172 A->factor = FACTOR_CHOLESKY; 173 PLogFlops((mat->n*mat->n*mat->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 (!mat->m || !mat->n) PetscFunctionReturn(0); 187 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 188 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 189 ierr = PetscMemcpy(y,x,mat->m*sizeof(Scalar));CHKERRQ(ierr); 190 if (A->factor == FACTOR_LU) { 191 LAgetrs_("N",&mat->m,&one,mat->v,&mat->m,mat->pivots,y,&mat->m,&info); 192 } else if (A->factor == FACTOR_CHOLESKY){ 193 LApotrs_("L",&mat->m,&one,mat->v,&mat->m,y,&mat->m,&info); 194 } 195 else SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Matrix must be factored to solve"); 196 if (info) SETERRQ(PETSC_ERR_LIB,0,"MBad solve"); 197 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 198 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 199 PLogFlops(2*mat->n*mat->n - mat->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 (!mat->m || !mat->n) PetscFunctionReturn(0); 213 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 214 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 215 ierr = PetscMemcpy(y,x,mat->m*sizeof(Scalar));CHKERRQ(ierr); 216 /* assume if pivots exist then use LU; else Cholesky */ 217 if (mat->pivots) { 218 LAgetrs_("T",&mat->m,&one,mat->v,&mat->m,mat->pivots,y,&mat->m,&info); 219 } else { 220 LApotrs_("L",&mat->m,&one,mat->v,&mat->m,y,&mat->m,&info); 221 } 222 if (info) SETERRQ(PETSC_ERR_LIB,0,"Bad solve"); 223 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 224 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 225 PLogFlops(2*mat->n*mat->n - mat->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 (!mat->m || !mat->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,mat->m*sizeof(Scalar));CHKERRQ(ierr); 248 /* assume if pivots exist then use LU; else Cholesky */ 249 if (mat->pivots) { 250 LAgetrs_("N",&mat->m,&one,mat->v,&mat->m,mat->pivots,y,&mat->m,&info); 251 } else { 252 LApotrs_("L",&mat->m,&one,mat->v,&mat->m,y,&mat->m,&info); 253 } 254 if (info) SETERRQ(PETSC_ERR_LIB,0,"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*mat->n*mat->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 (!mat->m || !mat->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,mat->m*sizeof(Scalar));CHKERRQ(ierr); 282 /* assume if pivots exist then use LU; else Cholesky */ 283 if (mat->pivots) { 284 LAgetrs_("T",&mat->m,&one,mat->v,&mat->m,mat->pivots,y,&mat->m,&info); 285 } else { 286 LApotrs_("L",&mat->m,&one,mat->v,&mat->m,y,&mat->m,&info); 287 } 288 if (info) SETERRQ(PETSC_ERR_LIB,0,"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*mat->n*mat->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 = mat->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 (!mat->m || !mat->n) PetscFunctionReturn(0); 372 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 373 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 374 LAgemv_("T",&(mat->m),&(mat->n),&_DOne,v,&(mat->m),x,&_One,&_DZero,y,&_One); 375 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 376 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 377 PLogFlops(2*mat->m*mat->n - mat->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 (!mat->m || !mat->n) PetscFunctionReturn(0); 391 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 392 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 393 LAgemv_("N",&(mat->m),&(mat->n),&_DOne,v,&(mat->m),x,&_One,&_DZero,y,&_One); 394 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 395 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 396 PLogFlops(2*mat->m*mat->n - mat->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 (!mat->m || !mat->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",&(mat->m),&(mat->n),&_DOne,v,&(mat->m),x,&_One,&_DOne,y,&_One); 414 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 415 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 416 PLogFlops(2*mat->m*mat->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 (!mat->m || !mat->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",&(mat->m),&(mat->n),&_DOne,v,&(mat->m),x,&_One,&_DOne,y,&_One); 435 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 436 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 437 PLogFlops(2*mat->m*mat->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 = mat->n; 452 if (cols) { 453 *cols = (int*)PetscMalloc((mat->n+1)*sizeof(int));CHKPTRQ(*cols); 454 for (i=0; i<mat->n; i++) (*cols)[i] = i; 455 } 456 if (vals) { 457 *vals = (Scalar*)PetscMalloc((mat->n+1)*sizeof(Scalar));CHKPTRQ(*vals); 458 v = mat->v + row; 459 for (i=0; i<mat->n; i++) {(*vals)[i] = *v; v += mat->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,0,"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,0,"Row too large"); 495 #endif 496 mat->v[indexn[j]*mat->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,0,"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,0,"Row too large"); 509 #endif 510 mat->v[indexn[j]*mat->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,0,"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,0,"Column too large"); 525 #endif 526 mat->v[indexn[j]*mat->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,0,"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,0,"Column too large"); 539 #endif 540 mat->v[indexn[j]*mat->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]*mat->m + indexm[i]]; 561 } 562 } 563 PetscFunctionReturn(0); 564 } 565 566 /* -----------------------------------------------------------------*/ 567 568 #include "sys.h" 569 570 #undef __FUNC__ 571 #define __FUNC__ /*<a name=""></a>*/"MatLoad_SeqDense" 572 int MatLoad_SeqDense(Viewer viewer,MatType type,Mat *A) 573 { 574 Mat_SeqDense *a; 575 Mat B; 576 int *scols,i,j,nz,ierr,fd,header[4],size; 577 int *rowlengths = 0,M,N,*cols; 578 Scalar *vals,*svals,*v,*w; 579 MPI_Comm comm = ((PetscObject)viewer)->comm; 580 581 PetscFunctionBegin; 582 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 583 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,0,"view must have one processor"); 584 ierr = ViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 585 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 586 if (header[0] != MAT_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,0,"Not matrix object"); 587 M = header[1]; N = header[2]; nz = header[3]; 588 589 if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense */ 590 ierr = MatCreateSeqDense(comm,M,N,PETSC_NULL,A);CHKERRQ(ierr); 591 B = *A; 592 a = (Mat_SeqDense*)B->data; 593 v = a->v; 594 /* Allocate some temp space to read in the values and then flip them 595 from row major to column major */ 596 w = (Scalar*)PetscMalloc((M*N+1)*sizeof(Scalar));CHKPTRQ(w); 597 /* read in nonzero values */ 598 ierr = PetscBinaryRead(fd,w,M*N,PETSC_SCALAR);CHKERRQ(ierr); 599 /* now flip the values and store them in the matrix*/ 600 for (j=0; j<N; j++) { 601 for (i=0; i<M; i++) { 602 *v++ =w[i*N+j]; 603 } 604 } 605 ierr = PetscFree(w);CHKERRQ(ierr); 606 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 607 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 608 } else { 609 /* read row lengths */ 610 rowlengths = (int*)PetscMalloc((M+1)*sizeof(int));CHKPTRQ(rowlengths); 611 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 612 613 /* create our matrix */ 614 ierr = MatCreateSeqDense(comm,M,N,PETSC_NULL,A);CHKERRQ(ierr); 615 B = *A; 616 a = (Mat_SeqDense*)B->data; 617 v = a->v; 618 619 /* read column indices and nonzeros */ 620 cols = scols = (int*)PetscMalloc((nz+1)*sizeof(int));CHKPTRQ(cols); 621 ierr = PetscBinaryRead(fd,cols,nz,PETSC_INT);CHKERRQ(ierr); 622 vals = svals = (Scalar*)PetscMalloc((nz+1)*sizeof(Scalar));CHKPTRQ(vals); 623 ierr = PetscBinaryRead(fd,vals,nz,PETSC_SCALAR);CHKERRQ(ierr); 624 625 /* insert into matrix */ 626 for (i=0; i<M; i++) { 627 for (j=0; j<rowlengths[i]; j++) v[i+M*scols[j]] = svals[j]; 628 svals += rowlengths[i]; scols += rowlengths[i]; 629 } 630 ierr = PetscFree(vals);CHKERRQ(ierr); 631 ierr = PetscFree(cols);CHKERRQ(ierr); 632 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 633 634 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 635 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 636 } 637 PetscFunctionReturn(0); 638 } 639 640 #include "sys.h" 641 642 #undef __FUNC__ 643 #define __FUNC__ /*<a name=""></a>*/"MatView_SeqDense_ASCII" 644 static int MatView_SeqDense_ASCII(Mat A,Viewer viewer) 645 { 646 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 647 int ierr,i,j,format; 648 char *outputname; 649 Scalar *v; 650 651 PetscFunctionBegin; 652 ierr = ViewerGetOutputname(viewer,&outputname);CHKERRQ(ierr); 653 ierr = ViewerGetFormat(viewer,&format);CHKERRQ(ierr); 654 if (format == VIEWER_FORMAT_ASCII_INFO || format == VIEWER_FORMAT_ASCII_INFO_LONG) { 655 PetscFunctionReturn(0); /* do nothing for now */ 656 } else if (format == VIEWER_FORMAT_ASCII_COMMON) { 657 ierr = ViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 658 for (i=0; i<a->m; i++) { 659 v = a->v + i; 660 ierr = ViewerASCIIPrintf(viewer,"row %d:",i);CHKERRQ(ierr); 661 for (j=0; j<a->n; j++) { 662 #if defined(PETSC_USE_COMPLEX) 663 if (PetscRealPart(*v) != 0.0 && PetscImaginaryPart(*v) != 0.0) { 664 ierr = ViewerASCIIPrintf(viewer," %d %g + %g i",j,PetscRealPart(*v),PetscImaginaryPart(*v));CHKERRQ(ierr); 665 } else if (PetscRealPart(*v)) { 666 ierr = ViewerASCIIPrintf(viewer," %d %g ",j,PetscRealPart(*v));CHKERRQ(ierr); 667 } 668 #else 669 if (*v) { 670 ierr = ViewerASCIIPrintf(viewer," %d %g ",j,*v);CHKERRQ(ierr); 671 } 672 #endif 673 v += a->m; 674 } 675 ierr = ViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 676 } 677 ierr = ViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 678 } else { 679 ierr = ViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 680 #if defined(PETSC_USE_COMPLEX) 681 int allreal = 1; 682 /* determine if matrix has all real values */ 683 v = a->v; 684 for (i=0; i<a->m*a->n; i++) { 685 if (PetscImaginaryPart(v[i])) { allreal = 0; break ;} 686 } 687 #endif 688 for (i=0; i<a->m; i++) { 689 v = a->v + i; 690 for (j=0; j<a->n; j++) { 691 #if defined(PETSC_USE_COMPLEX) 692 if (allreal) { 693 ierr = ViewerASCIIPrintf(viewer,"%6.4e ",PetscRealPart(*v));CHKERRQ(ierr); v += a->m; 694 } else { 695 ierr = ViewerASCIIPrintf(viewer,"%6.4e + %6.4e i ",PetscRealPart(*v),PetscImaginaryPart(*v));CHKERRQ(ierr); v += a->m; 696 } 697 #else 698 ierr = ViewerASCIIPrintf(viewer,"%6.4e ",*v);CHKERRQ(ierr); v += a->m; 699 #endif 700 } 701 ierr = ViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 702 } 703 ierr = ViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 704 } 705 ierr = ViewerFlush(viewer);CHKERRQ(ierr); 706 PetscFunctionReturn(0); 707 } 708 709 #undef __FUNC__ 710 #define __FUNC__ /*<a name=""></a>*/"MatView_SeqDense_Binary" 711 static int MatView_SeqDense_Binary(Mat A,Viewer viewer) 712 { 713 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 714 int ict,j,n = a->n,m = a->m,i,fd,*col_lens,ierr,nz = m*n; 715 int format; 716 Scalar *v,*anonz,*vals; 717 718 PetscFunctionBegin; 719 ierr = ViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 720 721 ierr = ViewerGetFormat(viewer,&format);CHKERRQ(ierr); 722 if (format == VIEWER_FORMAT_BINARY_NATIVE) { 723 /* store the matrix as a dense matrix */ 724 col_lens = (int*)PetscMalloc(4*sizeof(int));CHKPTRQ(col_lens); 725 col_lens[0] = MAT_COOKIE; 726 col_lens[1] = m; 727 col_lens[2] = n; 728 col_lens[3] = MATRIX_BINARY_FORMAT_DENSE; 729 ierr = PetscBinaryWrite(fd,col_lens,4,PETSC_INT,1);CHKERRQ(ierr); 730 ierr = PetscFree(col_lens);CHKERRQ(ierr); 731 732 /* write out matrix, by rows */ 733 vals = (Scalar*)PetscMalloc((m*n+1)*sizeof(Scalar));CHKPTRQ(vals); 734 v = a->v; 735 for (i=0; i<m; i++) { 736 for (j=0; j<n; j++) { 737 vals[i + j*m] = *v++; 738 } 739 } 740 ierr = PetscBinaryWrite(fd,vals,n*m,PETSC_SCALAR,0);CHKERRQ(ierr); 741 ierr = PetscFree(vals);CHKERRQ(ierr); 742 } else { 743 col_lens = (int*)PetscMalloc((4+nz)*sizeof(int));CHKPTRQ(col_lens); 744 col_lens[0] = MAT_COOKIE; 745 col_lens[1] = m; 746 col_lens[2] = n; 747 col_lens[3] = nz; 748 749 /* store lengths of each row and write (including header) to file */ 750 for (i=0; i<m; i++) col_lens[4+i] = n; 751 ierr = PetscBinaryWrite(fd,col_lens,4+m,PETSC_INT,1);CHKERRQ(ierr); 752 753 /* Possibly should write in smaller increments, not whole matrix at once? */ 754 /* store column indices (zero start index) */ 755 ict = 0; 756 for (i=0; i<m; i++) { 757 for (j=0; j<n; j++) col_lens[ict++] = j; 758 } 759 ierr = PetscBinaryWrite(fd,col_lens,nz,PETSC_INT,0);CHKERRQ(ierr); 760 ierr = PetscFree(col_lens);CHKERRQ(ierr); 761 762 /* store nonzero values */ 763 anonz = (Scalar*)PetscMalloc((nz+1)*sizeof(Scalar));CHKPTRQ(anonz); 764 ict = 0; 765 for (i=0; i<m; i++) { 766 v = a->v + i; 767 for (j=0; j<n; j++) { 768 anonz[ict++] = *v; v += a->m; 769 } 770 } 771 ierr = PetscBinaryWrite(fd,anonz,nz,PETSC_SCALAR,0);CHKERRQ(ierr); 772 ierr = PetscFree(anonz);CHKERRQ(ierr); 773 } 774 PetscFunctionReturn(0); 775 } 776 777 #undef __FUNC__ 778 #define __FUNC__ /*<a name=""></a>*/"MatView_SeqDense_Draw_Zoom" 779 int MatView_SeqDense_Draw_Zoom(Draw draw,void *Aa) 780 { 781 Mat A = (Mat) Aa; 782 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 783 int m = a->m,n = a->n,format,color,i,j,ierr; 784 Scalar *v = a->v; 785 Viewer viewer; 786 Draw popup; 787 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r,scale,maxv = 0.0; 788 789 PetscFunctionBegin; 790 791 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 792 ierr = ViewerGetFormat(viewer,&format);CHKERRQ(ierr); 793 ierr = DrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 794 795 /* Loop over matrix elements drawing boxes */ 796 if (format != VIEWER_FORMAT_DRAW_CONTOUR) { 797 /* Blue for negative and Red for positive */ 798 color = DRAW_BLUE; 799 for(j = 0; j < n; j++) { 800 x_l = j; 801 x_r = x_l + 1.0; 802 for(i = 0; i < m; i++) { 803 y_l = m - i - 1.0; 804 y_r = y_l + 1.0; 805 #if defined(PETSC_USE_COMPLEX) 806 if (PetscRealPart(v[j*m+i]) > 0.) { 807 color = DRAW_RED; 808 } else if (PetscRealPart(v[j*m+i]) < 0.) { 809 color = DRAW_BLUE; 810 } else { 811 continue; 812 } 813 #else 814 if (v[j*m+i] > 0.) { 815 color = DRAW_RED; 816 } else if (v[j*m+i] < 0.) { 817 color = DRAW_BLUE; 818 } else { 819 continue; 820 } 821 #endif 822 ierr = DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 823 } 824 } 825 } else { 826 /* use contour shading to indicate magnitude of values */ 827 /* first determine max of all nonzero values */ 828 for(i = 0; i < m*n; i++) { 829 if (PetscAbsScalar(v[i]) > maxv) maxv = PetscAbsScalar(v[i]); 830 } 831 scale = (245.0 - DRAW_BASIC_COLORS)/maxv; 832 ierr = DrawGetPopup(draw,&popup);CHKERRQ(ierr); 833 if (popup) {ierr = DrawScalePopup(popup,0.0,maxv);CHKERRQ(ierr);} 834 for(j = 0; j < n; j++) { 835 x_l = j; 836 x_r = x_l + 1.0; 837 for(i = 0; i < m; i++) { 838 y_l = m - i - 1.0; 839 y_r = y_l + 1.0; 840 color = DRAW_BASIC_COLORS + (int)(scale*PetscAbsScalar(v[j*m+i])); 841 ierr = DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 842 } 843 } 844 } 845 PetscFunctionReturn(0); 846 } 847 848 #undef __FUNC__ 849 #define __FUNC__ /*<a name=""></a>*/"MatView_SeqDense_Draw" 850 int MatView_SeqDense_Draw(Mat A,Viewer viewer) 851 { 852 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 853 Draw draw; 854 PetscTruth isnull; 855 PetscReal xr,yr,xl,yl,h,w; 856 int ierr; 857 858 PetscFunctionBegin; 859 ierr = ViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 860 ierr = DrawIsNull(draw,&isnull);CHKERRQ(ierr); 861 if (isnull == PETSC_TRUE) PetscFunctionReturn(0); 862 863 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 864 xr = a->n; yr = a->m; h = yr/10.0; w = xr/10.0; 865 xr += w; yr += h; xl = -w; yl = -h; 866 ierr = DrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 867 ierr = DrawZoom(draw,MatView_SeqDense_Draw_Zoom,A);CHKERRQ(ierr); 868 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);CHKERRQ(ierr); 869 PetscFunctionReturn(0); 870 } 871 872 #undef __FUNC__ 873 #define __FUNC__ /*<a name=""></a>*/"MatView_SeqDense" 874 int MatView_SeqDense(Mat A,Viewer viewer) 875 { 876 Mat_SeqDense *a = (Mat_SeqDense*)A->data; 877 int ierr; 878 PetscTruth issocket,isascii,isbinary,isdraw; 879 880 PetscFunctionBegin; 881 ierr = PetscTypeCompare((PetscObject)viewer,SOCKET_VIEWER,&issocket);CHKERRQ(ierr); 882 ierr = PetscTypeCompare((PetscObject)viewer,ASCII_VIEWER,&isascii);CHKERRQ(ierr); 883 ierr = PetscTypeCompare((PetscObject)viewer,BINARY_VIEWER,&isbinary);CHKERRQ(ierr); 884 ierr = PetscTypeCompare((PetscObject)viewer,DRAW_VIEWER,&isdraw);CHKERRQ(ierr); 885 886 if (issocket) { 887 ierr = ViewerSocketPutScalar_Private(viewer,a->m,a->n,a->v);CHKERRQ(ierr); 888 } else if (isascii) { 889 ierr = MatView_SeqDense_ASCII(A,viewer);CHKERRQ(ierr); 890 } else if (isbinary) { 891 ierr = MatView_SeqDense_Binary(A,viewer);CHKERRQ(ierr); 892 } else if (isdraw) { 893 ierr = MatView_SeqDense_Draw(A,viewer);CHKERRQ(ierr); 894 } else { 895 SETERRQ1(1,1,"Viewer type %s not supported by dense matrix",((PetscObject)viewer)->type_name); 896 } 897 PetscFunctionReturn(0); 898 } 899 900 #undef __FUNC__ 901 #define __FUNC__ /*<a name=""></a>*/"MatDestroy_SeqDense" 902 int MatDestroy_SeqDense(Mat mat) 903 { 904 Mat_SeqDense *l = (Mat_SeqDense*)mat->data; 905 int ierr; 906 907 PetscFunctionBegin; 908 909 if (mat->mapping) { 910 ierr = ISLocalToGlobalMappingDestroy(mat->mapping);CHKERRQ(ierr); 911 } 912 if (mat->bmapping) { 913 ierr = ISLocalToGlobalMappingDestroy(mat->bmapping);CHKERRQ(ierr); 914 } 915 #if defined(PETSC_USE_LOG) 916 PLogObjectState((PetscObject)mat,"Rows %d Cols %d",l->m,l->n); 917 #endif 918 if (l->pivots) {ierr = PetscFree(l->pivots);CHKERRQ(ierr);} 919 if (!l->user_alloc) {ierr = PetscFree(l->v);CHKERRQ(ierr);} 920 ierr = PetscFree(l);CHKERRQ(ierr); 921 if (mat->rmap) { 922 ierr = MapDestroy(mat->rmap);CHKERRQ(ierr); 923 } 924 if (mat->cmap) { 925 ierr = MapDestroy(mat->cmap);CHKERRQ(ierr); 926 } 927 PLogObjectDestroy(mat); 928 PetscHeaderDestroy(mat); 929 PetscFunctionReturn(0); 930 } 931 932 #undef __FUNC__ 933 #define __FUNC__ /*<a name=""></a>*/"MatTranspose_SeqDense" 934 int MatTranspose_SeqDense(Mat A,Mat *matout) 935 { 936 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 937 int k,j,m,n,ierr; 938 Scalar *v,tmp; 939 940 PetscFunctionBegin; 941 v = mat->v; m = mat->m; n = mat->n; 942 if (!matout) { /* in place transpose */ 943 if (m != n) { /* malloc temp to hold transpose */ 944 Scalar *w = (Scalar*)PetscMalloc((m+1)*(n+1)*sizeof(Scalar));CHKPTRQ(w); 945 for (j=0; j<m; j++) { 946 for (k=0; k<n; k++) { 947 w[k + j*n] = v[j + k*m]; 948 } 949 } 950 ierr = PetscMemcpy(v,w,m*n*sizeof(Scalar));CHKERRQ(ierr); 951 ierr = PetscFree(w);CHKERRQ(ierr); 952 } else { 953 for (j=0; j<m; j++) { 954 for (k=0; k<j; k++) { 955 tmp = v[j + k*n]; 956 v[j + k*n] = v[k + j*n]; 957 v[k + j*n] = tmp; 958 } 959 } 960 } 961 } else { /* out-of-place transpose */ 962 Mat tmat; 963 Mat_SeqDense *tmatd; 964 Scalar *v2; 965 ierr = MatCreateSeqDense(A->comm,mat->n,mat->m,PETSC_NULL,&tmat);CHKERRQ(ierr); 966 tmatd = (Mat_SeqDense*)tmat->data; 967 v = mat->v; v2 = tmatd->v; 968 for (j=0; j<n; j++) { 969 for (k=0; k<m; k++) v2[j + k*n] = v[k + j*m]; 970 } 971 ierr = MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 972 ierr = MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 973 *matout = tmat; 974 } 975 PetscFunctionReturn(0); 976 } 977 978 #undef __FUNC__ 979 #define __FUNC__ /*<a name=""></a>*/"MatEqual_SeqDense" 980 int MatEqual_SeqDense(Mat A1,Mat A2,PetscTruth *flg) 981 { 982 Mat_SeqDense *mat1 = (Mat_SeqDense*)A1->data; 983 Mat_SeqDense *mat2 = (Mat_SeqDense*)A2->data; 984 int i; 985 Scalar *v1 = mat1->v,*v2 = mat2->v; 986 987 PetscFunctionBegin; 988 if (A2->type != MATSEQDENSE) SETERRQ(PETSC_ERR_SUP,0,"Matrices should be of type MATSEQDENSE"); 989 if (mat1->m != mat2->m) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 990 if (mat1->n != mat2->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 991 for (i=0; i<mat1->m*mat1->n; i++) { 992 if (*v1 != *v2) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 993 v1++; v2++; 994 } 995 *flg = PETSC_TRUE; 996 PetscFunctionReturn(0); 997 } 998 999 #undef __FUNC__ 1000 #define __FUNC__ /*<a name=""></a>*/"MatGetDiagonal_SeqDense" 1001 int MatGetDiagonal_SeqDense(Mat A,Vec v) 1002 { 1003 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1004 int ierr,i,n,len; 1005 Scalar *x,zero = 0.0; 1006 1007 PetscFunctionBegin; 1008 ierr = VecSet(&zero,v);CHKERRQ(ierr); 1009 ierr = VecGetSize(v,&n);CHKERRQ(ierr); 1010 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1011 len = PetscMin(mat->m,mat->n); 1012 if (n != mat->m) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Nonconforming mat and vec"); 1013 for (i=0; i<len; i++) { 1014 x[i] = mat->v[i*mat->m + i]; 1015 } 1016 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1017 PetscFunctionReturn(0); 1018 } 1019 1020 #undef __FUNC__ 1021 #define __FUNC__ /*<a name=""></a>*/"MatDiagonalScale_SeqDense" 1022 int MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr) 1023 { 1024 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1025 Scalar *l,*r,x,*v; 1026 int ierr,i,j,m = mat->m,n = mat->n; 1027 1028 PetscFunctionBegin; 1029 if (ll) { 1030 ierr = VecGetSize(ll,&m);CHKERRQ(ierr); 1031 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1032 if (m != mat->m) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Left scaling vec wrong size"); 1033 for (i=0; i<m; i++) { 1034 x = l[i]; 1035 v = mat->v + i; 1036 for (j=0; j<n; j++) { (*v) *= x; v+= m;} 1037 } 1038 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1039 PLogFlops(n*m); 1040 } 1041 if (rr) { 1042 ierr = VecGetSize(rr,&n);CHKERRQ(ierr); 1043 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1044 if (n != mat->n) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Right scaling vec wrong size"); 1045 for (i=0; i<n; i++) { 1046 x = r[i]; 1047 v = mat->v + i*m; 1048 for (j=0; j<m; j++) { (*v++) *= x;} 1049 } 1050 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1051 PLogFlops(n*m); 1052 } 1053 PetscFunctionReturn(0); 1054 } 1055 1056 #undef __FUNC__ 1057 #define __FUNC__ /*<a name=""></a>*/"MatNorm_SeqDense" 1058 int MatNorm_SeqDense(Mat A,NormType type,PetscReal *norm) 1059 { 1060 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1061 Scalar *v = mat->v; 1062 PetscReal sum = 0.0; 1063 int i,j; 1064 1065 PetscFunctionBegin; 1066 if (type == NORM_FROBENIUS) { 1067 for (i=0; i<mat->n*mat->m; i++) { 1068 #if defined(PETSC_USE_COMPLEX) 1069 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1070 #else 1071 sum += (*v)*(*v); v++; 1072 #endif 1073 } 1074 *norm = sqrt(sum); 1075 PLogFlops(2*mat->n*mat->m); 1076 } else if (type == NORM_1) { 1077 *norm = 0.0; 1078 for (j=0; j<mat->n; j++) { 1079 sum = 0.0; 1080 for (i=0; i<mat->m; i++) { 1081 sum += PetscAbsScalar(*v); v++; 1082 } 1083 if (sum > *norm) *norm = sum; 1084 } 1085 PLogFlops(mat->n*mat->m); 1086 } else if (type == NORM_INFINITY) { 1087 *norm = 0.0; 1088 for (j=0; j<mat->m; j++) { 1089 v = mat->v + j; 1090 sum = 0.0; 1091 for (i=0; i<mat->n; i++) { 1092 sum += PetscAbsScalar(*v); v += mat->m; 1093 } 1094 if (sum > *norm) *norm = sum; 1095 } 1096 PLogFlops(mat->n*mat->m); 1097 } else { 1098 SETERRQ(PETSC_ERR_SUP,0,"No two norm"); 1099 } 1100 PetscFunctionReturn(0); 1101 } 1102 1103 #undef __FUNC__ 1104 #define __FUNC__ /*<a name=""></a>*/"MatSetOption_SeqDense" 1105 int MatSetOption_SeqDense(Mat A,MatOption op) 1106 { 1107 Mat_SeqDense *aij = (Mat_SeqDense*)A->data; 1108 1109 PetscFunctionBegin; 1110 if (op == MAT_ROW_ORIENTED) aij->roworiented = 1; 1111 else if (op == MAT_COLUMN_ORIENTED) aij->roworiented = 0; 1112 else if (op == MAT_ROWS_SORTED || 1113 op == MAT_ROWS_UNSORTED || 1114 op == MAT_COLUMNS_SORTED || 1115 op == MAT_COLUMNS_UNSORTED || 1116 op == MAT_SYMMETRIC || 1117 op == MAT_STRUCTURALLY_SYMMETRIC || 1118 op == MAT_NO_NEW_NONZERO_LOCATIONS || 1119 op == MAT_YES_NEW_NONZERO_LOCATIONS || 1120 op == MAT_NEW_NONZERO_LOCATION_ERR || 1121 op == MAT_NO_NEW_DIAGONALS || 1122 op == MAT_YES_NEW_DIAGONALS || 1123 op == MAT_IGNORE_OFF_PROC_ENTRIES || 1124 op == MAT_USE_HASH_TABLE) { 1125 PLogInfo(A,"MatSetOption_SeqDense:Option ignored\n"); 1126 } else { 1127 SETERRQ(PETSC_ERR_SUP,0,"unknown option"); 1128 } 1129 PetscFunctionReturn(0); 1130 } 1131 1132 #undef __FUNC__ 1133 #define __FUNC__ /*<a name=""></a>*/"MatZeroEntries_SeqDense" 1134 int MatZeroEntries_SeqDense(Mat A) 1135 { 1136 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1137 int ierr; 1138 1139 PetscFunctionBegin; 1140 ierr = PetscMemzero(l->v,l->m*l->n*sizeof(Scalar));CHKERRQ(ierr); 1141 PetscFunctionReturn(0); 1142 } 1143 1144 #undef __FUNC__ 1145 #define __FUNC__ /*<a name=""></a>*/"MatGetBlockSize_SeqDense" 1146 int MatGetBlockSize_SeqDense(Mat A,int *bs) 1147 { 1148 PetscFunctionBegin; 1149 *bs = 1; 1150 PetscFunctionReturn(0); 1151 } 1152 1153 #undef __FUNC__ 1154 #define __FUNC__ /*<a name=""></a>*/"MatZeroRows_SeqDense" 1155 int MatZeroRows_SeqDense(Mat A,IS is,Scalar *diag) 1156 { 1157 Mat_SeqDense *l = (Mat_SeqDense*)A->data; 1158 int n = l->n,i,j,ierr,N,*rows; 1159 Scalar *slot; 1160 1161 PetscFunctionBegin; 1162 ierr = ISGetSize(is,&N);CHKERRQ(ierr); 1163 ierr = ISGetIndices(is,&rows);CHKERRQ(ierr); 1164 for (i=0; i<N; i++) { 1165 slot = l->v + rows[i]; 1166 for (j=0; j<n; j++) { *slot = 0.0; slot += n;} 1167 } 1168 if (diag) { 1169 for (i=0; i<N; i++) { 1170 slot = l->v + (n+1)*rows[i]; 1171 *slot = *diag; 1172 } 1173 } 1174 ISRestoreIndices(is,&rows); 1175 PetscFunctionReturn(0); 1176 } 1177 1178 #undef __FUNC__ 1179 #define __FUNC__ /*<a name=""></a>*/"MatGetSize_SeqDense" 1180 int MatGetSize_SeqDense(Mat A,int *m,int *n) 1181 { 1182 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1183 1184 PetscFunctionBegin; 1185 if (m) *m = mat->m; 1186 if (n) *n = mat->n; 1187 PetscFunctionReturn(0); 1188 } 1189 1190 #undef __FUNC__ 1191 #define __FUNC__ /*<a name=""></a>*/"MatGetOwnershipRange_SeqDense" 1192 int MatGetOwnershipRange_SeqDense(Mat A,int *m,int *n) 1193 { 1194 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1195 1196 PetscFunctionBegin; 1197 if (m) *m = 0; 1198 if (n) *n = mat->m; 1199 PetscFunctionReturn(0); 1200 } 1201 1202 #undef __FUNC__ 1203 #define __FUNC__ /*<a name=""></a>*/"MatGetArray_SeqDense" 1204 int MatGetArray_SeqDense(Mat A,Scalar **array) 1205 { 1206 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1207 1208 PetscFunctionBegin; 1209 *array = mat->v; 1210 PetscFunctionReturn(0); 1211 } 1212 1213 #undef __FUNC__ 1214 #define __FUNC__ /*<a name=""></a>*/"MatRestoreArray_SeqDense" 1215 int MatRestoreArray_SeqDense(Mat A,Scalar **array) 1216 { 1217 PetscFunctionBegin; 1218 *array = 0; /* user cannot accidently use the array later */ 1219 PetscFunctionReturn(0); 1220 } 1221 1222 #undef __FUNC__ 1223 #define __FUNC__ /*<a name=""></a>*/"MatGetSubMatrix_SeqDense" 1224 static int MatGetSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,int cs,MatReuse scall,Mat *B) 1225 { 1226 Mat_SeqDense *mat = (Mat_SeqDense*)A->data; 1227 int i,j,ierr,m = mat->m,*irow,*icol,nrows,ncols; 1228 Scalar *av,*bv,*v = mat->v; 1229 Mat newmat; 1230 1231 PetscFunctionBegin; 1232 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1233 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1234 ierr = ISGetSize(isrow,&nrows);CHKERRQ(ierr); 1235 ierr = ISGetSize(iscol,&ncols);CHKERRQ(ierr); 1236 1237 /* Check submatrixcall */ 1238 if (scall == MAT_REUSE_MATRIX) { 1239 int n_cols,n_rows; 1240 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1241 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Reused submatrix wrong size"); 1242 newmat = *B; 1243 } else { 1244 /* Create and fill new matrix */ 1245 ierr = MatCreateSeqDense(A->comm,nrows,ncols,PETSC_NULL,&newmat);CHKERRQ(ierr); 1246 } 1247 1248 /* Now extract the data pointers and do the copy,column at a time */ 1249 bv = ((Mat_SeqDense*)newmat->data)->v; 1250 1251 for (i=0; i<ncols; i++) { 1252 av = v + m*icol[i]; 1253 for (j=0; j<nrows; j++) { 1254 *bv++ = av[irow[j]]; 1255 } 1256 } 1257 1258 /* Assemble the matrices so that the correct flags are set */ 1259 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1260 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1261 1262 /* Free work space */ 1263 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1264 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1265 *B = newmat; 1266 PetscFunctionReturn(0); 1267 } 1268 1269 #undef __FUNC__ 1270 #define __FUNC__ /*<a name=""></a>*/"MatGetSubMatrices_SeqDense" 1271 int MatGetSubMatrices_SeqDense(Mat A,int n,IS *irow,IS *icol,MatReuse scall,Mat **B) 1272 { 1273 int ierr,i; 1274 1275 PetscFunctionBegin; 1276 if (scall == MAT_INITIAL_MATRIX) { 1277 *B = (Mat*)PetscMalloc((n+1)*sizeof(Mat));CHKPTRQ(*B); 1278 } 1279 1280 for (i=0; i<n; i++) { 1281 ierr = MatGetSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1282 } 1283 PetscFunctionReturn(0); 1284 } 1285 1286 #undef __FUNC__ 1287 #define __FUNC__ /*<a name=""></a>*/"MatCopy_SeqDense" 1288 int MatCopy_SeqDense(Mat A,Mat B,MatStructure str) 1289 { 1290 Mat_SeqDense *a = (Mat_SeqDense*)A->data,*b = (Mat_SeqDense *)B->data; 1291 int ierr; 1292 1293 PetscFunctionBegin; 1294 if (B->type != MATSEQDENSE) { 1295 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1296 PetscFunctionReturn(0); 1297 } 1298 if (a->m != b->m || a->n != b->n) SETERRQ(PETSC_ERR_ARG_SIZ,0,"size(B) != size(A)"); 1299 ierr = PetscMemcpy(b->v,a->v,a->m*a->n*sizeof(Scalar));CHKERRQ(ierr); 1300 PetscFunctionReturn(0); 1301 } 1302 1303 /* -------------------------------------------------------------------*/ 1304 static struct _MatOps MatOps_Values = {MatSetValues_SeqDense, 1305 MatGetRow_SeqDense, 1306 MatRestoreRow_SeqDense, 1307 MatMult_SeqDense, 1308 MatMultAdd_SeqDense, 1309 MatMultTranspose_SeqDense, 1310 MatMultTransposeAdd_SeqDense, 1311 MatSolve_SeqDense, 1312 MatSolveAdd_SeqDense, 1313 MatSolveTranspose_SeqDense, 1314 MatSolveTransposeAdd_SeqDense, 1315 MatLUFactor_SeqDense, 1316 MatCholeskyFactor_SeqDense, 1317 MatRelax_SeqDense, 1318 MatTranspose_SeqDense, 1319 MatGetInfo_SeqDense, 1320 MatEqual_SeqDense, 1321 MatGetDiagonal_SeqDense, 1322 MatDiagonalScale_SeqDense, 1323 MatNorm_SeqDense, 1324 0, 1325 0, 1326 0, 1327 MatSetOption_SeqDense, 1328 MatZeroEntries_SeqDense, 1329 MatZeroRows_SeqDense, 1330 MatLUFactorSymbolic_SeqDense, 1331 MatLUFactorNumeric_SeqDense, 1332 MatCholeskyFactorSymbolic_SeqDense, 1333 MatCholeskyFactorNumeric_SeqDense, 1334 MatGetSize_SeqDense, 1335 MatGetSize_SeqDense, 1336 MatGetOwnershipRange_SeqDense, 1337 0, 1338 0, 1339 MatGetArray_SeqDense, 1340 MatRestoreArray_SeqDense, 1341 MatDuplicate_SeqDense, 1342 0, 1343 0, 1344 0, 1345 0, 1346 MatAXPY_SeqDense, 1347 MatGetSubMatrices_SeqDense, 1348 0, 1349 MatGetValues_SeqDense, 1350 MatCopy_SeqDense, 1351 0, 1352 MatScale_SeqDense, 1353 0, 1354 0, 1355 0, 1356 MatGetBlockSize_SeqDense, 1357 0, 1358 0, 1359 0, 1360 0, 1361 0, 1362 0, 1363 0, 1364 0, 1365 0, 1366 0, 1367 0, 1368 0, 1369 MatGetMaps_Petsc}; 1370 1371 #undef __FUNC__ 1372 #define __FUNC__ /*<a name=""></a>*/"MatCreateSeqDense" 1373 /*@C 1374 MatCreateSeqDense - Creates a sequential dense matrix that 1375 is stored in column major order (the usual Fortran 77 manner). Many 1376 of the matrix operations use the BLAS and LAPACK routines. 1377 1378 Collective on MPI_Comm 1379 1380 Input Parameters: 1381 + comm - MPI communicator, set to PETSC_COMM_SELF 1382 . m - number of rows 1383 . n - number of columns 1384 - data - optional location of matrix data. Set data=PETSC_NULL for PETSc 1385 to control all matrix memory allocation. 1386 1387 Output Parameter: 1388 . A - the matrix 1389 1390 Notes: 1391 The data input variable is intended primarily for Fortran programmers 1392 who wish to allocate their own matrix memory space. Most users should 1393 set data=PETSC_NULL. 1394 1395 Level: intermediate 1396 1397 .keywords: dense, matrix, LAPACK, BLAS 1398 1399 .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues() 1400 @*/ 1401 int MatCreateSeqDense(MPI_Comm comm,int m,int n,Scalar *data,Mat *A) 1402 { 1403 Mat B; 1404 Mat_SeqDense *b; 1405 int ierr,size; 1406 PetscTruth flg; 1407 1408 PetscFunctionBegin; 1409 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 1410 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,0,"Comm must be of size 1"); 1411 1412 *A = 0; 1413 PetscHeaderCreate(B,_p_Mat,struct _MatOps,MAT_COOKIE,MATSEQDENSE,"Mat",comm,MatDestroy,MatView); 1414 PLogObjectCreate(B); 1415 b = (Mat_SeqDense*)PetscMalloc(sizeof(Mat_SeqDense));CHKPTRQ(b); 1416 ierr = PetscMemzero(b,sizeof(Mat_SeqDense));CHKERRQ(ierr); 1417 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1418 B->ops->destroy = MatDestroy_SeqDense; 1419 B->ops->view = MatView_SeqDense; 1420 B->factor = 0; 1421 B->mapping = 0; 1422 PLogObjectMemory(B,sizeof(struct _p_Mat)); 1423 B->data = (void*)b; 1424 1425 b->m = m; B->m = m; B->M = m; 1426 b->n = n; B->n = n; B->N = n; 1427 1428 ierr = MapCreateMPI(comm,m,m,&B->rmap);CHKERRQ(ierr); 1429 ierr = MapCreateMPI(comm,n,n,&B->cmap);CHKERRQ(ierr); 1430 1431 b->pivots = 0; 1432 b->roworiented = 1; 1433 if (!data) { 1434 b->v = (Scalar*)PetscMalloc((m*n+1)*sizeof(Scalar));CHKPTRQ(b->v); 1435 ierr = PetscMemzero(b->v,m*n*sizeof(Scalar));CHKERRQ(ierr); 1436 b->user_alloc = 0; 1437 PLogObjectMemory(B,n*m*sizeof(Scalar)); 1438 } else { /* user-allocated storage */ 1439 b->v = data; 1440 b->user_alloc = 1; 1441 } 1442 ierr = OptionsHasName(PETSC_NULL,"-help",&flg);CHKERRQ(ierr); 1443 if (flg) {ierr = MatPrintHelp(B);CHKERRQ(ierr);} 1444 1445 *A = B; 1446 PetscFunctionReturn(0); 1447 } 1448 1449 1450 1451 1452 1453 1454