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