1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: dense.c,v 1.157 1998/10/01 18:54:20 bsmith Exp balay $"; 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__ "MatDuplicate_SeqDense" 94 int MatDuplicate_SeqDense(Mat A,MatDuplicateOption cpvalues,Mat *newmat) 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 == MAT_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 = MatDuplicate_SeqDense(A,MAT_DO_NOT_COPY_VALUES,fact);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,*w; 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 v = a->v; 598 /* Allocate some temp space to read in the values and then flip them 599 from row major to column major */ 600 w = (Scalar *) PetscMalloc((M*N+1)*sizeof(Scalar));CHKPTRQ(w); 601 /* read in nonzero values */ 602 ierr = PetscBinaryRead(fd,w,M*N,PETSC_SCALAR); CHKERRQ(ierr); 603 /* now flip the values and store them in the matrix*/ 604 for ( j=0; j<N; j++ ) { 605 for ( i=0; i<M; i++ ) { 606 *v++ =w[i*N+j]; 607 } 608 } 609 PetscFree(w); 610 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 611 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 612 } else { 613 /* read row lengths */ 614 rowlengths = (int*) PetscMalloc( (M+1)*sizeof(int) ); CHKPTRQ(rowlengths); 615 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT); CHKERRQ(ierr); 616 617 /* create our matrix */ 618 ierr = MatCreateSeqDense(comm,M,N,PETSC_NULL,A); CHKERRQ(ierr); 619 B = *A; 620 a = (Mat_SeqDense *) B->data; 621 v = a->v; 622 623 /* read column indices and nonzeros */ 624 cols = scols = (int *) PetscMalloc( (nz+1)*sizeof(int) ); CHKPTRQ(cols); 625 ierr = PetscBinaryRead(fd,cols,nz,PETSC_INT); CHKERRQ(ierr); 626 vals = svals = (Scalar *) PetscMalloc( (nz+1)*sizeof(Scalar) ); CHKPTRQ(vals); 627 ierr = PetscBinaryRead(fd,vals,nz,PETSC_SCALAR); CHKERRQ(ierr); 628 629 /* insert into matrix */ 630 for ( i=0; i<M; i++ ) { 631 for ( j=0; j<rowlengths[i]; j++ ) v[i+M*scols[j]] = svals[j]; 632 svals += rowlengths[i]; scols += rowlengths[i]; 633 } 634 PetscFree(vals); PetscFree(cols); PetscFree(rowlengths); 635 636 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 637 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 638 } 639 PetscFunctionReturn(0); 640 } 641 642 #include "pinclude/pviewer.h" 643 #include "sys.h" 644 645 #undef __FUNC__ 646 #define __FUNC__ "MatView_SeqDense_ASCII" 647 static int MatView_SeqDense_ASCII(Mat A,Viewer viewer) 648 { 649 Mat_SeqDense *a = (Mat_SeqDense *) A->data; 650 int ierr, i, j, format; 651 FILE *fd; 652 char *outputname; 653 Scalar *v; 654 655 PetscFunctionBegin; 656 ierr = ViewerASCIIGetPointer(viewer,&fd); CHKERRQ(ierr); 657 ierr = ViewerFileGetOutputname_Private(viewer,&outputname); CHKERRQ(ierr); 658 ierr = ViewerGetFormat(viewer,&format); 659 if (format == VIEWER_FORMAT_ASCII_INFO || format == VIEWER_FORMAT_ASCII_INFO_LONG) { 660 PetscFunctionReturn(0); /* do nothing for now */ 661 } 662 else if (format == VIEWER_FORMAT_ASCII_COMMON) { 663 for ( i=0; i<a->m; i++ ) { 664 v = a->v + i; 665 fprintf(fd,"row %d:",i); 666 for ( j=0; j<a->n; j++ ) { 667 #if defined(USE_PETSC_COMPLEX) 668 if (PetscReal(*v) != 0.0 && PetscImaginary(*v) != 0.0) fprintf(fd," %d %g + %g i",j,PetscReal(*v),PetscImaginary(*v)); 669 else if (PetscReal(*v)) fprintf(fd," %d %g ",j,PetscReal(*v)); 670 #else 671 if (*v) fprintf(fd," %d %g ",j,*v); 672 #endif 673 v += a->m; 674 } 675 fprintf(fd,"\n"); 676 } 677 } else { 678 #if defined(USE_PETSC_COMPLEX) 679 int allreal = 1; 680 /* determine if matrix has all real values */ 681 v = a->v; 682 for ( i=0; i<a->m*a->n; i++ ) { 683 if (PetscImaginary(v[i])) { allreal = 0; break ;} 684 } 685 #endif 686 for ( i=0; i<a->m; i++ ) { 687 v = a->v + i; 688 for ( j=0; j<a->n; j++ ) { 689 #if defined(USE_PETSC_COMPLEX) 690 if (allreal) { 691 fprintf(fd,"%6.4e ",PetscReal(*v)); v += a->m; 692 } else { 693 fprintf(fd,"%6.4e + %6.4e i ",PetscReal(*v),PetscImaginary(*v)); v += a->m; 694 } 695 #else 696 fprintf(fd,"%6.4e ",*v); v += a->m; 697 #endif 698 } 699 fprintf(fd,"\n"); 700 } 701 } 702 fflush(fd); 703 PetscFunctionReturn(0); 704 } 705 706 #undef __FUNC__ 707 #define __FUNC__ "MatView_SeqDense_Binary" 708 static int MatView_SeqDense_Binary(Mat A,Viewer viewer) 709 { 710 Mat_SeqDense *a = (Mat_SeqDense *) A->data; 711 int ict, j, n = a->n, m = a->m, i, fd, *col_lens, ierr, nz = m*n; 712 int format; 713 Scalar *v, *anonz,*vals; 714 715 PetscFunctionBegin; 716 ierr = ViewerBinaryGetDescriptor(viewer,&fd); CHKERRQ(ierr); 717 718 ierr = ViewerGetFormat(viewer,&format); CHKERRQ(ierr); 719 if (format == VIEWER_FORMAT_BINARY_NATIVE) { 720 /* store the matrix as a dense matrix */ 721 col_lens = (int *) PetscMalloc( 4*sizeof(int) ); CHKPTRQ(col_lens); 722 col_lens[0] = MAT_COOKIE; 723 col_lens[1] = m; 724 col_lens[2] = n; 725 col_lens[3] = MATRIX_BINARY_FORMAT_DENSE; 726 ierr = PetscBinaryWrite(fd,col_lens,4,PETSC_INT,1); CHKERRQ(ierr); 727 PetscFree(col_lens); 728 729 /* write out matrix, by rows */ 730 vals = (Scalar *) PetscMalloc((m*n+1)*sizeof(Scalar)); CHKPTRQ(vals); 731 v = a->v; 732 for ( i=0; i<m; i++ ) { 733 for ( j=0; j<n; j++ ) { 734 vals[i + j*m] = *v++; 735 } 736 } 737 ierr = PetscBinaryWrite(fd,vals,n*m,PETSC_SCALAR,0); CHKERRQ(ierr); 738 PetscFree(vals); 739 } else { 740 col_lens = (int *) PetscMalloc( (4+nz)*sizeof(int) ); CHKPTRQ(col_lens); 741 col_lens[0] = MAT_COOKIE; 742 col_lens[1] = m; 743 col_lens[2] = n; 744 col_lens[3] = nz; 745 746 /* store lengths of each row and write (including header) to file */ 747 for ( i=0; i<m; i++ ) col_lens[4+i] = n; 748 ierr = PetscBinaryWrite(fd,col_lens,4+m,PETSC_INT,1); CHKERRQ(ierr); 749 750 /* Possibly should write in smaller increments, not whole matrix at once? */ 751 /* store column indices (zero start index) */ 752 ict = 0; 753 for ( i=0; i<m; i++ ) { 754 for ( j=0; j<n; j++ ) col_lens[ict++] = j; 755 } 756 ierr = PetscBinaryWrite(fd,col_lens,nz,PETSC_INT,0); CHKERRQ(ierr); 757 PetscFree(col_lens); 758 759 /* store nonzero values */ 760 anonz = (Scalar *) PetscMalloc((nz+1)*sizeof(Scalar)); CHKPTRQ(anonz); 761 ict = 0; 762 for ( i=0; i<m; i++ ) { 763 v = a->v + i; 764 for ( j=0; j<n; j++ ) { 765 anonz[ict++] = *v; v += a->m; 766 } 767 } 768 ierr = PetscBinaryWrite(fd,anonz,nz,PETSC_SCALAR,0); CHKERRQ(ierr); 769 PetscFree(anonz); 770 } 771 PetscFunctionReturn(0); 772 } 773 774 #undef __FUNC__ 775 #define __FUNC__ "MatView_SeqDense" 776 int MatView_SeqDense(Mat A,Viewer viewer) 777 { 778 Mat_SeqDense *a = (Mat_SeqDense*) A->data; 779 ViewerType vtype; 780 int ierr; 781 782 PetscFunctionBegin; 783 ierr = ViewerGetType(viewer,&vtype); CHKERRQ(ierr); 784 785 if (vtype == MATLAB_VIEWER) { 786 ierr = ViewerMatlabPutScalar_Private(viewer,a->m,a->n,a->v); CHKERRQ(ierr); 787 } else if (vtype == ASCII_FILE_VIEWER || vtype == ASCII_FILES_VIEWER) { 788 ierr = MatView_SeqDense_ASCII(A,viewer);CHKERRQ(ierr); 789 } else if (vtype == BINARY_FILE_VIEWER) { 790 ierr = MatView_SeqDense_Binary(A,viewer);CHKERRQ(ierr); 791 } else { 792 SETERRQ(1,1,"Viewer type not supported by PETSc object"); 793 } 794 PetscFunctionReturn(0); 795 } 796 797 #undef __FUNC__ 798 #define __FUNC__ "MatDestroy_SeqDense" 799 int MatDestroy_SeqDense(Mat mat) 800 { 801 Mat_SeqDense *l = (Mat_SeqDense *) mat->data; 802 int ierr; 803 804 PetscFunctionBegin; 805 if (--mat->refct > 0) PetscFunctionReturn(0); 806 807 if (mat->mapping) { 808 ierr = ISLocalToGlobalMappingDestroy(mat->mapping); CHKERRQ(ierr); 809 } 810 if (mat->bmapping) { 811 ierr = ISLocalToGlobalMappingDestroy(mat->bmapping); CHKERRQ(ierr); 812 } 813 #if defined(USE_PETSC_LOG) 814 PLogObjectState((PetscObject)mat,"Rows %d Cols %d",l->m,l->n); 815 #endif 816 if (l->pivots) PetscFree(l->pivots); 817 if (!l->user_alloc) PetscFree(l->v); 818 PetscFree(l); 819 if (mat->rmap) { 820 ierr = MapDestroy(mat->rmap);CHKERRQ(ierr); 821 } 822 if (mat->cmap) { 823 ierr = MapDestroy(mat->cmap);CHKERRQ(ierr); 824 } 825 PLogObjectDestroy(mat); 826 PetscHeaderDestroy(mat); 827 PetscFunctionReturn(0); 828 } 829 830 #undef __FUNC__ 831 #define __FUNC__ "MatTranspose_SeqDense" 832 int MatTranspose_SeqDense(Mat A,Mat *matout) 833 { 834 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 835 int k, j, m, n; 836 Scalar *v, tmp; 837 838 PetscFunctionBegin; 839 v = mat->v; m = mat->m; n = mat->n; 840 if (matout == PETSC_NULL) { /* in place transpose */ 841 if (m != n) { /* malloc temp to hold transpose */ 842 Scalar *w = (Scalar *) PetscMalloc((m+1)*(n+1)*sizeof(Scalar));CHKPTRQ(w); 843 for ( j=0; j<m; j++ ) { 844 for ( k=0; k<n; k++ ) { 845 w[k + j*n] = v[j + k*m]; 846 } 847 } 848 PetscMemcpy(v,w,m*n*sizeof(Scalar)); 849 PetscFree(w); 850 } else { 851 for ( j=0; j<m; j++ ) { 852 for ( k=0; k<j; k++ ) { 853 tmp = v[j + k*n]; 854 v[j + k*n] = v[k + j*n]; 855 v[k + j*n] = tmp; 856 } 857 } 858 } 859 } else { /* out-of-place transpose */ 860 int ierr; 861 Mat tmat; 862 Mat_SeqDense *tmatd; 863 Scalar *v2; 864 ierr = MatCreateSeqDense(A->comm,mat->n,mat->m,PETSC_NULL,&tmat); CHKERRQ(ierr); 865 tmatd = (Mat_SeqDense *) tmat->data; 866 v = mat->v; v2 = tmatd->v; 867 for ( j=0; j<n; j++ ) { 868 for ( k=0; k<m; k++ ) v2[j + k*n] = v[k + j*m]; 869 } 870 ierr = MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 871 ierr = MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 872 *matout = tmat; 873 } 874 PetscFunctionReturn(0); 875 } 876 877 #undef __FUNC__ 878 #define __FUNC__ "MatEqual_SeqDense" 879 int MatEqual_SeqDense(Mat A1,Mat A2, PetscTruth *flg) 880 { 881 Mat_SeqDense *mat1 = (Mat_SeqDense *) A1->data; 882 Mat_SeqDense *mat2 = (Mat_SeqDense *) A2->data; 883 int i; 884 Scalar *v1 = mat1->v, *v2 = mat2->v; 885 886 PetscFunctionBegin; 887 if (A2->type != MATSEQDENSE) SETERRQ(PETSC_ERR_SUP,0,"Matrices should be of type MATSEQDENSE"); 888 if (mat1->m != mat2->m) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 889 if (mat1->n != mat2->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 890 for ( i=0; i<mat1->m*mat1->n; i++ ) { 891 if (*v1 != *v2) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 892 v1++; v2++; 893 } 894 *flg = PETSC_TRUE; 895 PetscFunctionReturn(0); 896 } 897 898 #undef __FUNC__ 899 #define __FUNC__ "MatGetDiagonal_SeqDense" 900 int MatGetDiagonal_SeqDense(Mat A,Vec v) 901 { 902 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 903 int ierr,i, n, len; 904 Scalar *x, zero = 0.0; 905 906 PetscFunctionBegin; 907 ierr = VecSet(&zero,v);CHKERRQ(ierr); 908 ierr = VecGetArray(v,&x); CHKERRQ(ierr); 909 ierr = VecGetSize(v,&n);CHKERRQ(ierr); 910 len = PetscMin(mat->m,mat->n); 911 if (n != mat->m) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Nonconforming mat and vec"); 912 for ( i=0; i<len; i++ ) { 913 x[i] = mat->v[i*mat->m + i]; 914 } 915 ierr = VecRestoreArray(v,&x); CHKERRQ(ierr); 916 PetscFunctionReturn(0); 917 } 918 919 #undef __FUNC__ 920 #define __FUNC__ "MatDiagonalScale_SeqDense" 921 int MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr) 922 { 923 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 924 Scalar *l,*r,x,*v; 925 int ierr,i,j,m = mat->m, n = mat->n; 926 927 PetscFunctionBegin; 928 if (ll) { 929 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 930 ierr = VecGetSize(ll,&m);CHKERRQ(ierr); 931 if (m != mat->m) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Left scaling vec wrong size"); 932 for ( i=0; i<m; i++ ) { 933 x = l[i]; 934 v = mat->v + i; 935 for ( j=0; j<n; j++ ) { (*v) *= x; v+= m;} 936 } 937 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 938 PLogFlops(n*m); 939 } 940 if (rr) { 941 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 942 ierr = VecGetSize(rr,&n);CHKERRQ(ierr); 943 if (n != mat->n) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Right scaling vec wrong size"); 944 for ( i=0; i<n; i++ ) { 945 x = r[i]; 946 v = mat->v + i*m; 947 for ( j=0; j<m; j++ ) { (*v++) *= x;} 948 } 949 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 950 PLogFlops(n*m); 951 } 952 PetscFunctionReturn(0); 953 } 954 955 #undef __FUNC__ 956 #define __FUNC__ "MatNorm_SeqDense" 957 int MatNorm_SeqDense(Mat A,NormType type,double *norm) 958 { 959 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 960 Scalar *v = mat->v; 961 double sum = 0.0; 962 int i, j; 963 964 PetscFunctionBegin; 965 if (type == NORM_FROBENIUS) { 966 for (i=0; i<mat->n*mat->m; i++ ) { 967 #if defined(USE_PETSC_COMPLEX) 968 sum += PetscReal(PetscConj(*v)*(*v)); v++; 969 #else 970 sum += (*v)*(*v); v++; 971 #endif 972 } 973 *norm = sqrt(sum); 974 PLogFlops(2*mat->n*mat->m); 975 } else if (type == NORM_1) { 976 *norm = 0.0; 977 for ( j=0; j<mat->n; j++ ) { 978 sum = 0.0; 979 for ( i=0; i<mat->m; i++ ) { 980 sum += PetscAbsScalar(*v); v++; 981 } 982 if (sum > *norm) *norm = sum; 983 } 984 PLogFlops(mat->n*mat->m); 985 } else if (type == NORM_INFINITY) { 986 *norm = 0.0; 987 for ( j=0; j<mat->m; j++ ) { 988 v = mat->v + j; 989 sum = 0.0; 990 for ( i=0; i<mat->n; i++ ) { 991 sum += PetscAbsScalar(*v); v += mat->m; 992 } 993 if (sum > *norm) *norm = sum; 994 } 995 PLogFlops(mat->n*mat->m); 996 } else { 997 SETERRQ(PETSC_ERR_SUP,0,"No two norm"); 998 } 999 PetscFunctionReturn(0); 1000 } 1001 1002 #undef __FUNC__ 1003 #define __FUNC__ "MatSetOption_SeqDense" 1004 int MatSetOption_SeqDense(Mat A,MatOption op) 1005 { 1006 Mat_SeqDense *aij = (Mat_SeqDense *) A->data; 1007 1008 PetscFunctionBegin; 1009 if (op == MAT_ROW_ORIENTED) aij->roworiented = 1; 1010 else if (op == MAT_COLUMN_ORIENTED) aij->roworiented = 0; 1011 else if (op == MAT_ROWS_SORTED || 1012 op == MAT_ROWS_UNSORTED || 1013 op == MAT_COLUMNS_SORTED || 1014 op == MAT_COLUMNS_UNSORTED || 1015 op == MAT_SYMMETRIC || 1016 op == MAT_STRUCTURALLY_SYMMETRIC || 1017 op == MAT_NO_NEW_NONZERO_LOCATIONS || 1018 op == MAT_YES_NEW_NONZERO_LOCATIONS || 1019 op == MAT_NEW_NONZERO_LOCATION_ERROR || 1020 op == MAT_NO_NEW_DIAGONALS || 1021 op == MAT_YES_NEW_DIAGONALS || 1022 op == MAT_IGNORE_OFF_PROC_ENTRIES || 1023 op == MAT_USE_HASH_TABLE) 1024 PLogInfo(A,"MatSetOption_SeqDense:Option ignored\n"); 1025 else { 1026 SETERRQ(PETSC_ERR_SUP,0,"unknown option"); 1027 } 1028 PetscFunctionReturn(0); 1029 } 1030 1031 #undef __FUNC__ 1032 #define __FUNC__ "MatZeroEntries_SeqDense" 1033 int MatZeroEntries_SeqDense(Mat A) 1034 { 1035 Mat_SeqDense *l = (Mat_SeqDense *) A->data; 1036 1037 PetscFunctionBegin; 1038 PetscMemzero(l->v,l->m*l->n*sizeof(Scalar)); 1039 PetscFunctionReturn(0); 1040 } 1041 1042 #undef __FUNC__ 1043 #define __FUNC__ "MatGetBlockSize_SeqDense" 1044 int MatGetBlockSize_SeqDense(Mat A,int *bs) 1045 { 1046 PetscFunctionBegin; 1047 *bs = 1; 1048 PetscFunctionReturn(0); 1049 } 1050 1051 #undef __FUNC__ 1052 #define __FUNC__ "MatZeroRows_SeqDense" 1053 int MatZeroRows_SeqDense(Mat A,IS is,Scalar *diag) 1054 { 1055 Mat_SeqDense *l = (Mat_SeqDense *) A->data; 1056 int n = l->n, i, j,ierr,N, *rows; 1057 Scalar *slot; 1058 1059 PetscFunctionBegin; 1060 ierr = ISGetSize(is,&N); CHKERRQ(ierr); 1061 ierr = ISGetIndices(is,&rows); CHKERRQ(ierr); 1062 for ( i=0; i<N; i++ ) { 1063 slot = l->v + rows[i]; 1064 for ( j=0; j<n; j++ ) { *slot = 0.0; slot += n;} 1065 } 1066 if (diag) { 1067 for ( i=0; i<N; i++ ) { 1068 slot = l->v + (n+1)*rows[i]; 1069 *slot = *diag; 1070 } 1071 } 1072 ISRestoreIndices(is,&rows); 1073 PetscFunctionReturn(0); 1074 } 1075 1076 #undef __FUNC__ 1077 #define __FUNC__ "MatGetSize_SeqDense" 1078 int MatGetSize_SeqDense(Mat A,int *m,int *n) 1079 { 1080 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 1081 1082 PetscFunctionBegin; 1083 *m = mat->m; *n = mat->n; 1084 PetscFunctionReturn(0); 1085 } 1086 1087 #undef __FUNC__ 1088 #define __FUNC__ "MatGetOwnershipRange_SeqDense" 1089 int MatGetOwnershipRange_SeqDense(Mat A,int *m,int *n) 1090 { 1091 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 1092 1093 PetscFunctionBegin; 1094 *m = 0; *n = mat->m; 1095 PetscFunctionReturn(0); 1096 } 1097 1098 #undef __FUNC__ 1099 #define __FUNC__ "MatGetArray_SeqDense" 1100 int MatGetArray_SeqDense(Mat A,Scalar **array) 1101 { 1102 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 1103 1104 PetscFunctionBegin; 1105 *array = mat->v; 1106 PetscFunctionReturn(0); 1107 } 1108 1109 #undef __FUNC__ 1110 #define __FUNC__ "MatRestoreArray_SeqDense" 1111 int MatRestoreArray_SeqDense(Mat A,Scalar **array) 1112 { 1113 PetscFunctionBegin; 1114 PetscFunctionReturn(0); 1115 } 1116 1117 #undef __FUNC__ 1118 #define __FUNC__ "MatGetSubMatrix_SeqDense" 1119 static int MatGetSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,int cs,MatGetSubMatrixCall scall,Mat *B) 1120 { 1121 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 1122 int i, j, ierr, m = mat->m, *irow, *icol, nrows, ncols; 1123 Scalar *av, *bv, *v = mat->v; 1124 Mat newmat; 1125 1126 PetscFunctionBegin; 1127 ierr = ISGetIndices(isrow,&irow); CHKERRQ(ierr); 1128 ierr = ISGetIndices(iscol,&icol); CHKERRQ(ierr); 1129 ierr = ISGetSize(isrow,&nrows); CHKERRQ(ierr); 1130 ierr = ISGetSize(iscol,&ncols); CHKERRQ(ierr); 1131 1132 /* Check submatrixcall */ 1133 if (scall == MAT_REUSE_MATRIX) { 1134 int n_cols,n_rows; 1135 ierr = MatGetSize(*B,&n_rows,&n_cols); CHKERRQ(ierr); 1136 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Reused submatrix wrong size"); 1137 newmat = *B; 1138 } else { 1139 /* Create and fill new matrix */ 1140 ierr = MatCreateSeqDense(A->comm,nrows,ncols,PETSC_NULL,&newmat); CHKERRQ(ierr); 1141 } 1142 1143 /* Now extract the data pointers and do the copy, column at a time */ 1144 bv = ((Mat_SeqDense*)newmat->data)->v; 1145 1146 for ( i=0; i<ncols; i++ ) { 1147 av = v + m*icol[i]; 1148 for (j=0; j<nrows; j++ ) { 1149 *bv++ = av[irow[j]]; 1150 } 1151 } 1152 1153 /* Assemble the matrices so that the correct flags are set */ 1154 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1155 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1156 1157 /* Free work space */ 1158 ierr = ISRestoreIndices(isrow,&irow); CHKERRQ(ierr); 1159 ierr = ISRestoreIndices(iscol,&icol); CHKERRQ(ierr); 1160 *B = newmat; 1161 PetscFunctionReturn(0); 1162 } 1163 1164 #undef __FUNC__ 1165 #define __FUNC__ "MatGetSubMatrices_SeqDense" 1166 int MatGetSubMatrices_SeqDense(Mat A,int n, IS *irow,IS *icol,MatGetSubMatrixCall scall,Mat **B) 1167 { 1168 int ierr,i; 1169 1170 PetscFunctionBegin; 1171 if (scall == MAT_INITIAL_MATRIX) { 1172 *B = (Mat *) PetscMalloc( (n+1)*sizeof(Mat) ); CHKPTRQ(*B); 1173 } 1174 1175 for ( i=0; i<n; i++ ) { 1176 ierr = MatGetSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1177 } 1178 PetscFunctionReturn(0); 1179 } 1180 1181 #undef __FUNC__ 1182 #define __FUNC__ "MatCopy_SeqDense" 1183 int MatCopy_SeqDense(Mat A, Mat B,MatStructure str) 1184 { 1185 Mat_SeqDense *a = (Mat_SeqDense *) A->data, *b = (Mat_SeqDense *)B->data; 1186 int ierr; 1187 1188 PetscFunctionBegin; 1189 if (B->type != MATSEQDENSE) { 1190 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1191 PetscFunctionReturn(0); 1192 } 1193 if (a->m != b->m || a->n != b->n) SETERRQ(PETSC_ERR_ARG_SIZ,0,"size(B) != size(A)"); 1194 PetscMemcpy(b->v,a->v,a->m*a->n*sizeof(Scalar)); 1195 PetscFunctionReturn(0); 1196 } 1197 1198 /* -------------------------------------------------------------------*/ 1199 static struct _MatOps MatOps_Values = {MatSetValues_SeqDense, 1200 MatGetRow_SeqDense, 1201 MatRestoreRow_SeqDense, 1202 MatMult_SeqDense, 1203 MatMultAdd_SeqDense, 1204 MatMultTrans_SeqDense, 1205 MatMultTransAdd_SeqDense, 1206 MatSolve_SeqDense, 1207 MatSolveAdd_SeqDense, 1208 MatSolveTrans_SeqDense, 1209 MatSolveTransAdd_SeqDense, 1210 MatLUFactor_SeqDense, 1211 MatCholeskyFactor_SeqDense, 1212 MatRelax_SeqDense, 1213 MatTranspose_SeqDense, 1214 MatGetInfo_SeqDense, 1215 MatEqual_SeqDense, 1216 MatGetDiagonal_SeqDense, 1217 MatDiagonalScale_SeqDense, 1218 MatNorm_SeqDense, 1219 0, 1220 0, 1221 0, 1222 MatSetOption_SeqDense, 1223 MatZeroEntries_SeqDense, 1224 MatZeroRows_SeqDense, 1225 MatLUFactorSymbolic_SeqDense, 1226 MatLUFactorNumeric_SeqDense, 1227 MatCholeskyFactorSymbolic_SeqDense, 1228 MatCholeskyFactorNumeric_SeqDense, 1229 MatGetSize_SeqDense, 1230 MatGetSize_SeqDense, 1231 MatGetOwnershipRange_SeqDense, 1232 0, 1233 0, 1234 MatGetArray_SeqDense, 1235 MatRestoreArray_SeqDense, 1236 MatDuplicate_SeqDense, 1237 0, 1238 0, 1239 0, 1240 0, 1241 MatAXPY_SeqDense, 1242 MatGetSubMatrices_SeqDense, 1243 0, 1244 MatGetValues_SeqDense, 1245 MatCopy_SeqDense, 1246 0, 1247 MatScale_SeqDense, 1248 0, 1249 0, 1250 0, 1251 MatGetBlockSize_SeqDense, 1252 0, 1253 0, 1254 0, 1255 0, 1256 0, 1257 0, 1258 0, 1259 0, 1260 0, 1261 0, 1262 0, 1263 0, 1264 MatGetMaps_Petsc}; 1265 1266 #undef __FUNC__ 1267 #define __FUNC__ "MatCreateSeqDense" 1268 /*@C 1269 MatCreateSeqDense - Creates a sequential dense matrix that 1270 is stored in column major order (the usual Fortran 77 manner). Many 1271 of the matrix operations use the BLAS and LAPACK routines. 1272 1273 Collective on MPI_Comm 1274 1275 Input Parameters: 1276 + comm - MPI communicator, set to PETSC_COMM_SELF 1277 . m - number of rows 1278 . n - number of columns 1279 - data - optional location of matrix data. Set data=PETSC_NULL for PETSc 1280 to control all matrix memory allocation. 1281 1282 Output Parameter: 1283 . A - the matrix 1284 1285 Notes: 1286 The data input variable is intended primarily for Fortran programmers 1287 who wish to allocate their own matrix memory space. Most users should 1288 set data=PETSC_NULL. 1289 1290 .keywords: dense, matrix, LAPACK, BLAS 1291 1292 .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues() 1293 @*/ 1294 int MatCreateSeqDense(MPI_Comm comm,int m,int n,Scalar *data,Mat *A) 1295 { 1296 Mat B; 1297 Mat_SeqDense *b; 1298 int ierr,flg,size; 1299 1300 PetscFunctionBegin; 1301 MPI_Comm_size(comm,&size); 1302 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,0,"Comm must be of size 1"); 1303 1304 *A = 0; 1305 PetscHeaderCreate(B,_p_Mat,struct _MatOps,MAT_COOKIE,MATSEQDENSE,comm,MatDestroy,MatView); 1306 PLogObjectCreate(B); 1307 b = (Mat_SeqDense *) PetscMalloc(sizeof(Mat_SeqDense)); CHKPTRQ(b); 1308 PetscMemzero(b,sizeof(Mat_SeqDense)); 1309 PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps)); 1310 B->ops->destroy = MatDestroy_SeqDense; 1311 B->ops->view = MatView_SeqDense; 1312 B->factor = 0; 1313 B->mapping = 0; 1314 PLogObjectMemory(B,sizeof(struct _p_Mat)); 1315 B->data = (void *) b; 1316 1317 b->m = m; B->m = m; B->M = m; 1318 b->n = n; B->n = n; B->N = n; 1319 1320 ierr = MapCreateMPI(comm,m,m,&B->rmap);CHKERRQ(ierr); 1321 ierr = MapCreateMPI(comm,n,n,&B->cmap);CHKERRQ(ierr); 1322 1323 b->pivots = 0; 1324 b->roworiented = 1; 1325 if (data == PETSC_NULL) { 1326 b->v = (Scalar*) PetscMalloc((m*n+1)*sizeof(Scalar)); CHKPTRQ(b->v); 1327 PetscMemzero(b->v,m*n*sizeof(Scalar)); 1328 b->user_alloc = 0; 1329 PLogObjectMemory(B,n*m*sizeof(Scalar)); 1330 } else { /* user-allocated storage */ 1331 b->v = data; 1332 b->user_alloc = 1; 1333 } 1334 ierr = OptionsHasName(PETSC_NULL,"-help",&flg); CHKERRQ(ierr); 1335 if (flg) {ierr = MatPrintHelp(B); CHKERRQ(ierr);} 1336 1337 *A = B; 1338 PetscFunctionReturn(0); 1339 } 1340 1341 1342 1343 1344 1345 1346