1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: dense.c,v 1.150 1998/05/29 22:43:06 balay 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__ "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->mapping) { 812 ierr = ISLocalToGlobalMappingDestroy(mat->mapping); CHKERRQ(ierr); 813 } 814 PLogObjectDestroy(mat); 815 PetscHeaderDestroy(mat); 816 PetscFunctionReturn(0); 817 } 818 819 #undef __FUNC__ 820 #define __FUNC__ "MatTranspose_SeqDense" 821 int MatTranspose_SeqDense(Mat A,Mat *matout) 822 { 823 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 824 int k, j, m, n; 825 Scalar *v, tmp; 826 827 PetscFunctionBegin; 828 v = mat->v; m = mat->m; n = mat->n; 829 if (matout == PETSC_NULL) { /* in place transpose */ 830 if (m != n) { /* malloc temp to hold transpose */ 831 Scalar *w = (Scalar *) PetscMalloc((m+1)*(n+1)*sizeof(Scalar));CHKPTRQ(w); 832 for ( j=0; j<m; j++ ) { 833 for ( k=0; k<n; k++ ) { 834 w[k + j*n] = v[j + k*m]; 835 } 836 } 837 PetscMemcpy(v,w,m*n*sizeof(Scalar)); 838 PetscFree(w); 839 } else { 840 for ( j=0; j<m; j++ ) { 841 for ( k=0; k<j; k++ ) { 842 tmp = v[j + k*n]; 843 v[j + k*n] = v[k + j*n]; 844 v[k + j*n] = tmp; 845 } 846 } 847 } 848 } else { /* out-of-place transpose */ 849 int ierr; 850 Mat tmat; 851 Mat_SeqDense *tmatd; 852 Scalar *v2; 853 ierr = MatCreateSeqDense(A->comm,mat->n,mat->m,PETSC_NULL,&tmat); CHKERRQ(ierr); 854 tmatd = (Mat_SeqDense *) tmat->data; 855 v = mat->v; v2 = tmatd->v; 856 for ( j=0; j<n; j++ ) { 857 for ( k=0; k<m; k++ ) v2[j + k*n] = v[k + j*m]; 858 } 859 ierr = MatAssemblyBegin(tmat,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 860 ierr = MatAssemblyEnd(tmat,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 861 *matout = tmat; 862 } 863 PetscFunctionReturn(0); 864 } 865 866 #undef __FUNC__ 867 #define __FUNC__ "MatEqual_SeqDense" 868 int MatEqual_SeqDense(Mat A1,Mat A2, PetscTruth *flg) 869 { 870 Mat_SeqDense *mat1 = (Mat_SeqDense *) A1->data; 871 Mat_SeqDense *mat2 = (Mat_SeqDense *) A2->data; 872 int i; 873 Scalar *v1 = mat1->v, *v2 = mat2->v; 874 875 PetscFunctionBegin; 876 if (A2->type != MATSEQDENSE) SETERRQ(PETSC_ERR_SUP,0,"Matrices should be of type MATSEQDENSE"); 877 if (mat1->m != mat2->m) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 878 if (mat1->n != mat2->n) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 879 for ( i=0; i<mat1->m*mat1->n; i++ ) { 880 if (*v1 != *v2) {*flg = PETSC_FALSE; PetscFunctionReturn(0);} 881 v1++; v2++; 882 } 883 *flg = PETSC_TRUE; 884 PetscFunctionReturn(0); 885 } 886 887 #undef __FUNC__ 888 #define __FUNC__ "MatGetDiagonal_SeqDense" 889 int MatGetDiagonal_SeqDense(Mat A,Vec v) 890 { 891 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 892 int ierr,i, n, len; 893 Scalar *x, zero = 0.0; 894 895 PetscFunctionBegin; 896 ierr = VecSet(&zero,v);CHKERRQ(ierr); 897 ierr = VecGetArray(v,&x); CHKERRQ(ierr); 898 ierr = VecGetSize(v,&n);CHKERRQ(ierr); 899 len = PetscMin(mat->m,mat->n); 900 if (n != mat->m) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Nonconforming mat and vec"); 901 for ( i=0; i<len; i++ ) { 902 x[i] = mat->v[i*mat->m + i]; 903 } 904 ierr = VecRestoreArray(v,&x); CHKERRQ(ierr); 905 PetscFunctionReturn(0); 906 } 907 908 #undef __FUNC__ 909 #define __FUNC__ "MatDiagonalScale_SeqDense" 910 int MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr) 911 { 912 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 913 Scalar *l,*r,x,*v; 914 int ierr,i,j,m = mat->m, n = mat->n; 915 916 PetscFunctionBegin; 917 if (ll) { 918 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 919 ierr = VecGetSize(ll,&m);CHKERRQ(ierr); 920 if (m != mat->m) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Left scaling vec wrong size"); 921 for ( i=0; i<m; i++ ) { 922 x = l[i]; 923 v = mat->v + i; 924 for ( j=0; j<n; j++ ) { (*v) *= x; v+= m;} 925 } 926 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 927 PLogFlops(n*m); 928 } 929 if (rr) { 930 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 931 ierr = VecGetSize(rr,&n);CHKERRQ(ierr); 932 if (n != mat->n) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Right scaling vec wrong size"); 933 for ( i=0; i<n; i++ ) { 934 x = r[i]; 935 v = mat->v + i*m; 936 for ( j=0; j<m; j++ ) { (*v++) *= x;} 937 } 938 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 939 PLogFlops(n*m); 940 } 941 PetscFunctionReturn(0); 942 } 943 944 #undef __FUNC__ 945 #define __FUNC__ "MatNorm_SeqDense" 946 int MatNorm_SeqDense(Mat A,NormType type,double *norm) 947 { 948 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 949 Scalar *v = mat->v; 950 double sum = 0.0; 951 int i, j; 952 953 PetscFunctionBegin; 954 if (type == NORM_FROBENIUS) { 955 for (i=0; i<mat->n*mat->m; i++ ) { 956 #if defined(USE_PETSC_COMPLEX) 957 sum += PetscReal(PetscConj(*v)*(*v)); v++; 958 #else 959 sum += (*v)*(*v); v++; 960 #endif 961 } 962 *norm = sqrt(sum); 963 PLogFlops(2*mat->n*mat->m); 964 } else if (type == NORM_1) { 965 *norm = 0.0; 966 for ( j=0; j<mat->n; j++ ) { 967 sum = 0.0; 968 for ( i=0; i<mat->m; i++ ) { 969 sum += PetscAbsScalar(*v); v++; 970 } 971 if (sum > *norm) *norm = sum; 972 } 973 PLogFlops(mat->n*mat->m); 974 } else if (type == NORM_INFINITY) { 975 *norm = 0.0; 976 for ( j=0; j<mat->m; j++ ) { 977 v = mat->v + j; 978 sum = 0.0; 979 for ( i=0; i<mat->n; i++ ) { 980 sum += PetscAbsScalar(*v); v += mat->m; 981 } 982 if (sum > *norm) *norm = sum; 983 } 984 PLogFlops(mat->n*mat->m); 985 } else { 986 SETERRQ(PETSC_ERR_SUP,0,"No two norm"); 987 } 988 PetscFunctionReturn(0); 989 } 990 991 #undef __FUNC__ 992 #define __FUNC__ "MatSetOption_SeqDense" 993 int MatSetOption_SeqDense(Mat A,MatOption op) 994 { 995 Mat_SeqDense *aij = (Mat_SeqDense *) A->data; 996 997 PetscFunctionBegin; 998 if (op == MAT_ROW_ORIENTED) aij->roworiented = 1; 999 else if (op == MAT_COLUMN_ORIENTED) aij->roworiented = 0; 1000 else if (op == MAT_ROWS_SORTED || 1001 op == MAT_ROWS_UNSORTED || 1002 op == MAT_COLUMNS_SORTED || 1003 op == MAT_COLUMNS_UNSORTED || 1004 op == MAT_SYMMETRIC || 1005 op == MAT_STRUCTURALLY_SYMMETRIC || 1006 op == MAT_NO_NEW_NONZERO_LOCATIONS || 1007 op == MAT_YES_NEW_NONZERO_LOCATIONS || 1008 op == MAT_NEW_NONZERO_LOCATION_ERROR || 1009 op == MAT_NO_NEW_DIAGONALS || 1010 op == MAT_YES_NEW_DIAGONALS || 1011 op == MAT_IGNORE_OFF_PROC_ENTRIES || 1012 op == MAT_USE_HASH_TABLE) 1013 PLogInfo(A,"MatSetOption_SeqDense:Option ignored\n"); 1014 else { 1015 SETERRQ(PETSC_ERR_SUP,0,"unknown option"); 1016 } 1017 PetscFunctionReturn(0); 1018 } 1019 1020 #undef __FUNC__ 1021 #define __FUNC__ "MatZeroEntries_SeqDense" 1022 int MatZeroEntries_SeqDense(Mat A) 1023 { 1024 Mat_SeqDense *l = (Mat_SeqDense *) A->data; 1025 1026 PetscFunctionBegin; 1027 PetscMemzero(l->v,l->m*l->n*sizeof(Scalar)); 1028 PetscFunctionReturn(0); 1029 } 1030 1031 #undef __FUNC__ 1032 #define __FUNC__ "MatGetBlockSize_SeqDense" 1033 int MatGetBlockSize_SeqDense(Mat A,int *bs) 1034 { 1035 PetscFunctionBegin; 1036 *bs = 1; 1037 PetscFunctionReturn(0); 1038 } 1039 1040 #undef __FUNC__ 1041 #define __FUNC__ "MatZeroRows_SeqDense" 1042 int MatZeroRows_SeqDense(Mat A,IS is,Scalar *diag) 1043 { 1044 Mat_SeqDense *l = (Mat_SeqDense *) A->data; 1045 int n = l->n, i, j,ierr,N, *rows; 1046 Scalar *slot; 1047 1048 PetscFunctionBegin; 1049 ierr = ISGetSize(is,&N); CHKERRQ(ierr); 1050 ierr = ISGetIndices(is,&rows); CHKERRQ(ierr); 1051 for ( i=0; i<N; i++ ) { 1052 slot = l->v + rows[i]; 1053 for ( j=0; j<n; j++ ) { *slot = 0.0; slot += n;} 1054 } 1055 if (diag) { 1056 for ( i=0; i<N; i++ ) { 1057 slot = l->v + (n+1)*rows[i]; 1058 *slot = *diag; 1059 } 1060 } 1061 ISRestoreIndices(is,&rows); 1062 PetscFunctionReturn(0); 1063 } 1064 1065 #undef __FUNC__ 1066 #define __FUNC__ "MatGetSize_SeqDense" 1067 int MatGetSize_SeqDense(Mat A,int *m,int *n) 1068 { 1069 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 1070 1071 PetscFunctionBegin; 1072 *m = mat->m; *n = mat->n; 1073 PetscFunctionReturn(0); 1074 } 1075 1076 #undef __FUNC__ 1077 #define __FUNC__ "MatGetOwnershipRange_SeqDense" 1078 int MatGetOwnershipRange_SeqDense(Mat A,int *m,int *n) 1079 { 1080 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 1081 1082 PetscFunctionBegin; 1083 *m = 0; *n = mat->m; 1084 PetscFunctionReturn(0); 1085 } 1086 1087 #undef __FUNC__ 1088 #define __FUNC__ "MatGetArray_SeqDense" 1089 int MatGetArray_SeqDense(Mat A,Scalar **array) 1090 { 1091 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 1092 1093 PetscFunctionBegin; 1094 *array = mat->v; 1095 PetscFunctionReturn(0); 1096 } 1097 1098 #undef __FUNC__ 1099 #define __FUNC__ "MatRestoreArray_SeqDense" 1100 int MatRestoreArray_SeqDense(Mat A,Scalar **array) 1101 { 1102 PetscFunctionBegin; 1103 PetscFunctionReturn(0); 1104 } 1105 1106 #undef __FUNC__ 1107 #define __FUNC__ "MatGetSubMatrix_SeqDense" 1108 static int MatGetSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,int cs,MatGetSubMatrixCall scall,Mat *submat) 1109 { 1110 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 1111 int nznew, *smap, i, j, ierr, oldcols = mat->n; 1112 int *irow, *icol, nrows, ncols, *cwork; 1113 Scalar *vwork, *val; 1114 Mat newmat; 1115 1116 PetscFunctionBegin; 1117 ierr = ISGetIndices(isrow,&irow); CHKERRQ(ierr); 1118 ierr = ISGetIndices(iscol,&icol); CHKERRQ(ierr); 1119 ierr = ISGetSize(isrow,&nrows); CHKERRQ(ierr); 1120 ierr = ISGetSize(iscol,&ncols); CHKERRQ(ierr); 1121 1122 smap = (int *) PetscMalloc(oldcols*sizeof(int)); CHKPTRQ(smap); 1123 cwork = (int *) PetscMalloc(ncols*sizeof(int)); CHKPTRQ(cwork); 1124 vwork = (Scalar *) PetscMalloc(ncols*sizeof(Scalar)); CHKPTRQ(vwork); 1125 PetscMemzero((char*)smap,oldcols*sizeof(int)); 1126 for ( i=0; i<ncols; i++ ) smap[icol[i]] = i+1; 1127 1128 /* Create and fill new matrix */ 1129 ierr = MatCreateSeqDense(A->comm,nrows,ncols,PETSC_NULL,&newmat); CHKERRQ(ierr); 1130 for (i=0; i<nrows; i++) { 1131 nznew = 0; 1132 val = mat->v + irow[i]; 1133 for (j=0; j<oldcols; j++) { 1134 if (smap[j]) { 1135 cwork[nznew] = smap[j] - 1; 1136 vwork[nznew++] = val[j * mat->m]; 1137 } 1138 } 1139 ierr = MatSetValues(newmat,1,&i,nznew,cwork,vwork,INSERT_VALUES);CHKERRQ(ierr); 1140 } 1141 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1142 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1143 1144 /* Free work space */ 1145 PetscFree(smap); PetscFree(cwork); PetscFree(vwork); 1146 ierr = ISRestoreIndices(isrow,&irow); CHKERRQ(ierr); 1147 ierr = ISRestoreIndices(iscol,&icol); CHKERRQ(ierr); 1148 *submat = newmat; 1149 PetscFunctionReturn(0); 1150 } 1151 1152 #undef __FUNC__ 1153 #define __FUNC__ "MatGetSubMatrices_SeqDense" 1154 int MatGetSubMatrices_SeqDense(Mat A,int n, IS *irow,IS *icol,MatGetSubMatrixCall scall,Mat **B) 1155 { 1156 int ierr,i; 1157 1158 PetscFunctionBegin; 1159 if (scall == MAT_INITIAL_MATRIX) { 1160 *B = (Mat *) PetscMalloc( (n+1)*sizeof(Mat) ); CHKPTRQ(*B); 1161 } 1162 1163 for ( i=0; i<n; i++ ) { 1164 ierr = MatGetSubMatrix_SeqDense(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1165 } 1166 PetscFunctionReturn(0); 1167 } 1168 1169 #undef __FUNC__ 1170 #define __FUNC__ "MatCopy_SeqDense" 1171 int MatCopy_SeqDense(Mat A, Mat B) 1172 { 1173 Mat_SeqDense *a = (Mat_SeqDense *) A->data, *b = (Mat_SeqDense *)B->data; 1174 int ierr; 1175 1176 PetscFunctionBegin; 1177 if (B->type != MATSEQDENSE) { 1178 ierr = MatCopy_Basic(A,B);CHKERRQ(ierr); 1179 PetscFunctionReturn(0); 1180 } 1181 if (a->m != b->m || a->n != b->n) SETERRQ(PETSC_ERR_ARG_SIZ,0,"size(B) != size(A)"); 1182 PetscMemcpy(b->v,a->v,a->m*a->n*sizeof(Scalar)); 1183 PetscFunctionReturn(0); 1184 } 1185 1186 /* -------------------------------------------------------------------*/ 1187 static struct _MatOps MatOps = {MatSetValues_SeqDense, 1188 MatGetRow_SeqDense, 1189 MatRestoreRow_SeqDense, 1190 MatMult_SeqDense, 1191 MatMultAdd_SeqDense, 1192 MatMultTrans_SeqDense, 1193 MatMultTransAdd_SeqDense, 1194 MatSolve_SeqDense, 1195 MatSolveAdd_SeqDense, 1196 MatSolveTrans_SeqDense, 1197 MatSolveTransAdd_SeqDense, 1198 MatLUFactor_SeqDense, 1199 MatCholeskyFactor_SeqDense, 1200 MatRelax_SeqDense, 1201 MatTranspose_SeqDense, 1202 MatGetInfo_SeqDense, 1203 MatEqual_SeqDense, 1204 MatGetDiagonal_SeqDense, 1205 MatDiagonalScale_SeqDense, 1206 MatNorm_SeqDense, 1207 0, 1208 0, 1209 0, 1210 MatSetOption_SeqDense, 1211 MatZeroEntries_SeqDense, 1212 MatZeroRows_SeqDense, 1213 MatLUFactorSymbolic_SeqDense, 1214 MatLUFactorNumeric_SeqDense, 1215 MatCholeskyFactorSymbolic_SeqDense, 1216 MatCholeskyFactorNumeric_SeqDense, 1217 MatGetSize_SeqDense, 1218 MatGetSize_SeqDense, 1219 MatGetOwnershipRange_SeqDense, 1220 0, 1221 0, 1222 MatGetArray_SeqDense, 1223 MatRestoreArray_SeqDense, 1224 MatConvertSameType_SeqDense,0,0,0,0, 1225 MatAXPY_SeqDense,MatGetSubMatrices_SeqDense,0, 1226 MatGetValues_SeqDense, 1227 MatCopy_SeqDense,0,MatScale_SeqDense, 1228 0,0,0,MatGetBlockSize_SeqDense}; 1229 1230 #undef __FUNC__ 1231 #define __FUNC__ "MatCreateSeqDense" 1232 /*@C 1233 MatCreateSeqDense - Creates a sequential dense matrix that 1234 is stored in column major order (the usual Fortran 77 manner). Many 1235 of the matrix operations use the BLAS and LAPACK routines. 1236 1237 Collective on MPI_Comm 1238 1239 Input Parameters: 1240 + comm - MPI communicator, set to PETSC_COMM_SELF 1241 . m - number of rows 1242 . n - number of columns 1243 - data - optional location of matrix data. Set data=PETSC_NULL for PETSc 1244 to control all matrix memory allocation. 1245 1246 Output Parameter: 1247 . A - the matrix 1248 1249 Notes: 1250 The data input variable is intended primarily for Fortran programmers 1251 who wish to allocate their own matrix memory space. Most users should 1252 set data=PETSC_NULL. 1253 1254 .keywords: dense, matrix, LAPACK, BLAS 1255 1256 .seealso: MatCreate(), MatCreateMPIDense(), MatSetValues() 1257 @*/ 1258 int MatCreateSeqDense(MPI_Comm comm,int m,int n,Scalar *data,Mat *A) 1259 { 1260 Mat B; 1261 Mat_SeqDense *b; 1262 int ierr,flg,size; 1263 1264 PetscFunctionBegin; 1265 MPI_Comm_size(comm,&size); 1266 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,0,"Comm must be of size 1"); 1267 1268 *A = 0; 1269 PetscHeaderCreate(B,_p_Mat,struct _MatOps,MAT_COOKIE,MATSEQDENSE,comm,MatDestroy,MatView); 1270 PLogObjectCreate(B); 1271 b = (Mat_SeqDense *) PetscMalloc(sizeof(Mat_SeqDense)); CHKPTRQ(b); 1272 PetscMemzero(b,sizeof(Mat_SeqDense)); 1273 PetscMemcpy(B->ops,&MatOps,sizeof(struct _MatOps)); 1274 B->ops->destroy = MatDestroy_SeqDense; 1275 B->ops->view = MatView_SeqDense; 1276 B->factor = 0; 1277 B->mapping = 0; 1278 PLogObjectMemory(B,sizeof(struct _p_Mat)); 1279 B->data = (void *) b; 1280 1281 b->m = m; B->m = m; B->M = m; 1282 b->n = n; B->n = n; B->N = n; 1283 b->pivots = 0; 1284 b->roworiented = 1; 1285 if (data == PETSC_NULL) { 1286 b->v = (Scalar*) PetscMalloc((m*n+1)*sizeof(Scalar)); CHKPTRQ(b->v); 1287 PetscMemzero(b->v,m*n*sizeof(Scalar)); 1288 b->user_alloc = 0; 1289 PLogObjectMemory(B,n*m*sizeof(Scalar)); 1290 } 1291 else { /* user-allocated storage */ 1292 b->v = data; 1293 b->user_alloc = 1; 1294 } 1295 ierr = OptionsHasName(PETSC_NULL,"-help",&flg); CHKERRQ(ierr); 1296 if (flg) {ierr = MatPrintHelp(B); CHKERRQ(ierr);} 1297 1298 *A = B; 1299 PetscFunctionReturn(0); 1300 } 1301 1302 1303 1304 1305 1306 1307