1 #ifndef lint 2 static char vcid[] = "$Id: dense.c,v 1.103 1996/04/09 02:23:22 curfman Exp balay $"; 3 #endif 4 /* 5 Defines the basic matrix operations for sequential dense. 6 */ 7 8 #include "dense.h" 9 #include "pinclude/plapack.h" 10 #include "pinclude/pviewer.h" 11 12 int MatAXPY_SeqDense(Scalar *alpha,Mat X,Mat Y) 13 { 14 Mat_SeqDense *x = (Mat_SeqDense*) X->data,*y = (Mat_SeqDense*) Y->data; 15 int N = x->m*x->n, one = 1; 16 BLaxpy_( &N, alpha, x->v, &one, y->v, &one ); 17 PLogFlops(2*N-1); 18 return 0; 19 } 20 21 static int MatGetInfo_SeqDense(Mat A,MatInfoType flag,int *nz,int *nzalloc,int *mem) 22 { 23 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 24 int i,N = mat->m*mat->n,count = 0; 25 Scalar *v = mat->v; 26 for ( i=0; i<N; i++ ) {if (*v != 0.0) count++; v++;} 27 if (nz) *nz = count; 28 if (nzalloc) *nzalloc = N; 29 if (mem) *mem = (int)A->mem; 30 return 0; 31 } 32 33 static int MatScale_SeqDense(Scalar *alpha,Mat inA) 34 { 35 Mat_SeqDense *a = (Mat_SeqDense *) inA->data; 36 int one = 1, nz; 37 38 nz = a->m*a->n; 39 BLscal_( &nz, alpha, a->v, &one ); 40 PLogFlops(nz); 41 return 0; 42 } 43 44 /* ---------------------------------------------------------------*/ 45 /* COMMENT: I have chosen to hide column permutation in the pivots, 46 rather than put it in the Mat->col slot.*/ 47 static int MatLUFactor_SeqDense(Mat A,IS row,IS col,double f) 48 { 49 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 50 int info; 51 52 if (!mat->pivots) { 53 mat->pivots = (int *) PetscMalloc(mat->m*sizeof(int));CHKPTRQ(mat->pivots); 54 PLogObjectMemory(A,mat->m*sizeof(int)); 55 } 56 LAgetrf_(&mat->m,&mat->n,mat->v,&mat->m,mat->pivots,&info); 57 if (info<0) SETERRQ(1,"MatLUFactor_SeqDense:Bad argument to LU factorization"); 58 if (info>0) SETERRQ(1,"MatLUFactor_SeqDense:Bad LU factorization"); 59 A->factor = FACTOR_LU; 60 PLogFlops((2*mat->n*mat->n*mat->n)/3); 61 return 0; 62 } 63 static int MatLUFactorSymbolic_SeqDense(Mat A,IS row,IS col,double f,Mat *fact) 64 { 65 return MatConvert(A,MATSAME,fact); 66 } 67 static int MatLUFactorNumeric_SeqDense(Mat A,Mat *fact) 68 { 69 return MatLUFactor(*fact,0,0,1.0); 70 } 71 static int MatCholeskyFactorSymbolic_SeqDense(Mat A,IS row,double f,Mat *fact) 72 { 73 return MatConvert(A,MATSAME,fact); 74 } 75 static int MatCholeskyFactorNumeric_SeqDense(Mat A,Mat *fact) 76 { 77 return MatCholeskyFactor(*fact,0,1.0); 78 } 79 static int MatCholeskyFactor_SeqDense(Mat A,IS perm,double f) 80 { 81 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 82 int info; 83 84 if (mat->pivots) { 85 PetscFree(mat->pivots); 86 PLogObjectMemory(A,-mat->m*sizeof(int)); 87 mat->pivots = 0; 88 } 89 LApotrf_("L",&mat->n,mat->v,&mat->m,&info); 90 if (info) SETERRQ(1,"MatCholeskyFactor_SeqDense:Bad factorization"); 91 A->factor = FACTOR_CHOLESKY; 92 PLogFlops((mat->n*mat->n*mat->n)/3); 93 return 0; 94 } 95 96 static int MatSolve_SeqDense(Mat A,Vec xx,Vec yy) 97 { 98 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 99 int one = 1, info, ierr; 100 Scalar *x, *y; 101 102 ierr = VecGetArray(xx,&x); CHKERRQ(ierr); ierr = VecGetArray(yy,&y); CHKERRQ(ierr); 103 PetscMemcpy(y,x,mat->m*sizeof(Scalar)); 104 if (A->factor == FACTOR_LU) { 105 LAgetrs_( "N", &mat->m, &one, mat->v, &mat->m, mat->pivots,y, &mat->m, &info ); 106 } 107 else if (A->factor == FACTOR_CHOLESKY){ 108 LApotrs_( "L", &mat->m, &one, mat->v, &mat->m,y, &mat->m, &info ); 109 } 110 else SETERRQ(1,"MatSolve_SeqDense:Matrix must be factored to solve"); 111 if (info) SETERRQ(1,"MatSolve_SeqDense:Bad solve"); 112 PLogFlops(mat->n*mat->n - mat->n); 113 return 0; 114 } 115 static int MatSolveTrans_SeqDense(Mat A,Vec xx,Vec yy) 116 { 117 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 118 int one = 1, info; 119 Scalar *x, *y; 120 121 VecGetArray(xx,&x); VecGetArray(yy,&y); 122 PetscMemcpy(y,x,mat->m*sizeof(Scalar)); 123 /* assume if pivots exist then use LU; else Cholesky */ 124 if (mat->pivots) { 125 LAgetrs_( "T", &mat->m, &one, mat->v, &mat->m, mat->pivots,y, &mat->m, &info ); 126 } 127 else { 128 LApotrs_( "L", &mat->m, &one, mat->v, &mat->m,y, &mat->m, &info ); 129 } 130 if (info) SETERRQ(1,"MatSolveTrans_SeqDense:Bad solve"); 131 PLogFlops(mat->n*mat->n - mat->n); 132 return 0; 133 } 134 static int MatSolveAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 135 { 136 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 137 int one = 1, info,ierr; 138 Scalar *x, *y, sone = 1.0; 139 Vec tmp = 0; 140 141 VecGetArray(xx,&x); VecGetArray(yy,&y); 142 if (yy == zz) { 143 ierr = VecDuplicate(yy,&tmp); CHKERRQ(ierr); 144 PLogObjectParent(A,tmp); 145 ierr = VecCopy(yy,tmp); CHKERRQ(ierr); 146 } 147 PetscMemcpy(y,x,mat->m*sizeof(Scalar)); 148 /* assume if pivots exist then use LU; else Cholesky */ 149 if (mat->pivots) { 150 LAgetrs_( "N", &mat->m, &one, mat->v, &mat->m, mat->pivots,y, &mat->m, &info ); 151 } 152 else { 153 LApotrs_( "L", &mat->m, &one, mat->v, &mat->m,y, &mat->m, &info ); 154 } 155 if (info) SETERRQ(1,"MatSolveAdd_SeqDense:Bad solve"); 156 if (tmp) {VecAXPY(&sone,tmp,yy); VecDestroy(tmp);} 157 else VecAXPY(&sone,zz,yy); 158 PLogFlops(mat->n*mat->n - mat->n); 159 return 0; 160 } 161 162 static int MatSolveTransAdd_SeqDense(Mat A,Vec xx,Vec zz, Vec yy) 163 { 164 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 165 int one = 1, info,ierr; 166 Scalar *x, *y, sone = 1.0; 167 Vec tmp; 168 169 VecGetArray(xx,&x); VecGetArray(yy,&y); 170 if (yy == zz) { 171 ierr = VecDuplicate(yy,&tmp); CHKERRQ(ierr); 172 PLogObjectParent(A,tmp); 173 ierr = VecCopy(yy,tmp); CHKERRQ(ierr); 174 } 175 PetscMemcpy(y,x,mat->m*sizeof(Scalar)); 176 /* assume if pivots exist then use LU; else Cholesky */ 177 if (mat->pivots) { 178 LAgetrs_( "T", &mat->m, &one, mat->v, &mat->m, mat->pivots,y, &mat->m, &info ); 179 } 180 else { 181 LApotrs_( "L", &mat->m, &one, mat->v, &mat->m,y, &mat->m, &info ); 182 } 183 if (info) SETERRQ(1,"MatSolveTransAdd_SeqDense:Bad solve"); 184 if (tmp) {VecAXPY(&sone,tmp,yy); VecDestroy(tmp);} 185 else VecAXPY(&sone,zz,yy); 186 PLogFlops(mat->n*mat->n - mat->n); 187 return 0; 188 } 189 /* ------------------------------------------------------------------*/ 190 static int MatRelax_SeqDense(Mat A,Vec bb,double omega,MatSORType flag, 191 double shift,int its,Vec xx) 192 { 193 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 194 Scalar *x, *b, *v = mat->v, zero = 0.0, xt; 195 int o = 1,ierr, m = mat->m, i; 196 197 if (flag & SOR_ZERO_INITIAL_GUESS) { 198 /* this is a hack fix, should have another version without 199 the second BLdot */ 200 ierr = VecSet(&zero,xx); CHKERRQ(ierr); 201 } 202 VecGetArray(xx,&x); VecGetArray(bb,&b); 203 while (its--) { 204 if (flag & SOR_FORWARD_SWEEP){ 205 for ( i=0; i<m; i++ ) { 206 #if defined(PETSC_COMPLEX) 207 /* cannot use BLAS dot for complex because compiler/linker is 208 not happy about returning a double complex */ 209 int _i; 210 Scalar sum = b[i]; 211 for ( _i=0; _i<m; _i++ ) { 212 sum -= conj(v[i+_i*m])*x[_i]; 213 } 214 xt = sum; 215 #else 216 xt = b[i]-BLdot_(&m,v+i,&m,x,&o); 217 #endif 218 x[i] = (1. - omega)*x[i] + omega*(xt/(v[i + i*m]+shift) + x[i]); 219 } 220 } 221 if (flag & SOR_BACKWARD_SWEEP) { 222 for ( i=m-1; i>=0; i-- ) { 223 #if defined(PETSC_COMPLEX) 224 /* cannot use BLAS dot for complex because compiler/linker is 225 not happy about returning a double complex */ 226 int _i; 227 Scalar sum = b[i]; 228 for ( _i=0; _i<m; _i++ ) { 229 sum -= conj(v[i+_i*m])*x[_i]; 230 } 231 xt = sum; 232 #else 233 xt = b[i]-BLdot_(&m,v+i,&m,x,&o); 234 #endif 235 x[i] = (1. - omega)*x[i] + omega*(xt/(v[i + i*m]+shift) + x[i]); 236 } 237 } 238 } 239 return 0; 240 } 241 242 /* -----------------------------------------------------------------*/ 243 int MatMultTrans_SeqDense(Mat A,Vec xx,Vec yy) 244 { 245 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 246 Scalar *v = mat->v, *x, *y; 247 int _One=1;Scalar _DOne=1.0, _DZero=0.0; 248 VecGetArray(xx,&x), VecGetArray(yy,&y); 249 LAgemv_("T",&(mat->m),&(mat->n),&_DOne,v,&(mat->m),x,&_One,&_DZero,y,&_One); 250 PLogFlops(2*mat->m*mat->n - mat->n); 251 return 0; 252 } 253 int MatMult_SeqDense(Mat A,Vec xx,Vec yy) 254 { 255 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 256 Scalar *v = mat->v, *x, *y; 257 int _One=1;Scalar _DOne=1.0, _DZero=0.0; 258 VecGetArray(xx,&x); VecGetArray(yy,&y); 259 LAgemv_( "N", &(mat->m), &(mat->n), &_DOne, v, &(mat->m),x,&_One,&_DZero,y,&_One); 260 PLogFlops(2*mat->m*mat->n - mat->m); 261 return 0; 262 } 263 int MatMultAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 264 { 265 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 266 Scalar *v = mat->v, *x, *y, *z; 267 int _One=1; Scalar _DOne=1.0; 268 VecGetArray(xx,&x); VecGetArray(yy,&y); VecGetArray(zz,&z); 269 if (zz != yy) PetscMemcpy(y,z,mat->m*sizeof(Scalar)); 270 LAgemv_( "N", &(mat->m), &(mat->n),&_DOne,v,&(mat->m),x,&_One,&_DOne,y,&_One); 271 PLogFlops(2*mat->m*mat->n); 272 return 0; 273 } 274 int MatMultTransAdd_SeqDense(Mat A,Vec xx,Vec zz,Vec yy) 275 { 276 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 277 Scalar *v = mat->v, *x, *y, *z; 278 int _One=1; 279 Scalar _DOne=1.0; 280 VecGetArray(xx,&x); VecGetArray(yy,&y); VecGetArray(zz,&z); 281 if (zz != yy) PetscMemcpy(y,z,mat->n*sizeof(Scalar)); 282 LAgemv_( "T", &(mat->m), &(mat->n), &_DOne, v, &(mat->m),x,&_One,&_DOne,y,&_One); 283 PLogFlops(2*mat->m*mat->n); 284 return 0; 285 } 286 287 /* -----------------------------------------------------------------*/ 288 static int MatGetRow_SeqDense(Mat A,int row,int *ncols,int **cols,Scalar **vals) 289 { 290 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 291 Scalar *v; 292 int i; 293 294 *ncols = mat->n; 295 if (cols) { 296 *cols = (int *) PetscMalloc(mat->n*sizeof(int)); CHKPTRQ(*cols); 297 for ( i=0; i<mat->n; i++ ) (*cols)[i] = i; 298 } 299 if (vals) { 300 *vals = (Scalar *) PetscMalloc(mat->n*sizeof(Scalar)); CHKPTRQ(*vals); 301 v = mat->v + row; 302 for ( i=0; i<mat->n; i++ ) {(*vals)[i] = *v; v += mat->m;} 303 } 304 return 0; 305 } 306 static int MatRestoreRow_SeqDense(Mat A,int row,int *ncols,int **cols,Scalar **vals) 307 { 308 if (cols) { PetscFree(*cols); } 309 if (vals) { PetscFree(*vals); } 310 return 0; 311 } 312 /* ----------------------------------------------------------------*/ 313 static int MatSetValues_SeqDense(Mat A,int m,int *indexm,int n, 314 int *indexn,Scalar *v,InsertMode addv) 315 { 316 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 317 int i,j; 318 319 if (!mat->roworiented) { 320 if (addv == INSERT_VALUES) { 321 for ( j=0; j<n; j++ ) { 322 for ( i=0; i<m; i++ ) { 323 mat->v[indexn[j]*mat->m + indexm[i]] = *v++; 324 } 325 } 326 } 327 else { 328 for ( j=0; j<n; j++ ) { 329 for ( i=0; i<m; i++ ) { 330 mat->v[indexn[j]*mat->m + indexm[i]] += *v++; 331 } 332 } 333 } 334 } 335 else { 336 if (addv == INSERT_VALUES) { 337 for ( i=0; i<m; i++ ) { 338 for ( j=0; j<n; j++ ) { 339 mat->v[indexn[j]*mat->m + indexm[i]] = *v++; 340 } 341 } 342 } 343 else { 344 for ( i=0; i<m; i++ ) { 345 for ( j=0; j<n; j++ ) { 346 mat->v[indexn[j]*mat->m + indexm[i]] += *v++; 347 } 348 } 349 } 350 } 351 return 0; 352 } 353 354 static int MatGetValues_SeqDense(Mat A,int m,int *indexm,int n,int *indexn,Scalar *v) 355 { 356 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 357 int i, j; 358 Scalar *vpt = v; 359 360 /* row-oriented output */ 361 for ( i=0; i<m; i++ ) { 362 for ( j=0; j<n; j++ ) { 363 *vpt++ = mat->v[indexn[j]*mat->m + indexm[i]]; 364 } 365 } 366 return 0; 367 } 368 369 /* -----------------------------------------------------------------*/ 370 static int MatConvertSameType_SeqDense(Mat A,Mat *newmat,int cpvalues) 371 { 372 Mat_SeqDense *mat = (Mat_SeqDense *) A->data, *l; 373 int ierr; 374 Mat newi; 375 376 ierr = MatCreateSeqDense(A->comm,mat->m,mat->n,PETSC_NULL,&newi); CHKERRQ(ierr); 377 l = (Mat_SeqDense *) newi->data; 378 if (cpvalues == COPY_VALUES) { 379 PetscMemcpy(l->v,mat->v,mat->m*mat->n*sizeof(Scalar)); 380 } 381 newi->assembled = PETSC_TRUE; 382 *newmat = newi; 383 return 0; 384 } 385 386 #include "sys.h" 387 388 int MatLoad_SeqDense(Viewer viewer,MatType type,Mat *A) 389 { 390 Mat_SeqDense *a; 391 Mat B; 392 int *scols, i, j, nz, ierr, fd, header[4], size; 393 int *rowlengths = 0, M, N, *cols; 394 Scalar *vals, *svals, *v; 395 MPI_Comm comm = ((PetscObject)viewer)->comm; 396 397 MPI_Comm_size(comm,&size); 398 if (size > 1) SETERRQ(1,"MatLoad_SeqDense: view must have one processor"); 399 ierr = ViewerBinaryGetDescriptor(viewer,&fd); CHKERRQ(ierr); 400 ierr = PetscBinaryRead(fd,header,4,BINARY_INT); CHKERRQ(ierr); 401 if (header[0] != MAT_COOKIE) SETERRQ(1,"MatLoad_SeqDense:Not matrix object"); 402 M = header[1]; N = header[2]; nz = header[3]; 403 404 if (nz == MATRIX_BINARY_FORMAT_DENSE) { /* matrix in file is dense */ 405 ierr = MatCreateSeqDense(comm,M,N,PETSC_NULL,A); CHKERRQ(ierr); 406 B = *A; 407 a = (Mat_SeqDense *) B->data; 408 409 /* read in nonzero values */ 410 ierr = PetscBinaryRead(fd,a->v,M*N,BINARY_SCALAR); CHKERRQ(ierr); 411 412 ierr = MatAssemblyBegin(B,FINAL_ASSEMBLY); CHKERRQ(ierr); 413 ierr = MatAssemblyEnd(B,FINAL_ASSEMBLY); CHKERRQ(ierr); 414 ierr = MatTranspose(B,PETSC_NULL); CHKERRQ(ierr); 415 } else { 416 /* read row lengths */ 417 rowlengths = (int*) PetscMalloc( M*sizeof(int) ); CHKPTRQ(rowlengths); 418 ierr = PetscBinaryRead(fd,rowlengths,M,BINARY_INT); CHKERRQ(ierr); 419 420 /* create our matrix */ 421 ierr = MatCreateSeqDense(comm,M,N,PETSC_NULL,A); CHKERRQ(ierr); 422 B = *A; 423 a = (Mat_SeqDense *) B->data; 424 v = a->v; 425 426 /* read column indices and nonzeros */ 427 cols = scols = (int *) PetscMalloc( nz*sizeof(int) ); CHKPTRQ(cols); 428 ierr = PetscBinaryRead(fd,cols,nz,BINARY_INT); CHKERRQ(ierr); 429 vals = svals = (Scalar *) PetscMalloc( nz*sizeof(Scalar) ); CHKPTRQ(vals); 430 ierr = PetscBinaryRead(fd,vals,nz,BINARY_SCALAR); CHKERRQ(ierr); 431 432 /* insert into matrix */ 433 for ( i=0; i<M; i++ ) { 434 for ( j=0; j<rowlengths[i]; j++ ) v[i+M*scols[j]] = svals[j]; 435 svals += rowlengths[i]; scols += rowlengths[i]; 436 } 437 PetscFree(vals); PetscFree(cols); PetscFree(rowlengths); 438 439 ierr = MatAssemblyBegin(B,FINAL_ASSEMBLY); CHKERRQ(ierr); 440 ierr = MatAssemblyEnd(B,FINAL_ASSEMBLY); CHKERRQ(ierr); 441 } 442 return 0; 443 } 444 445 #include "pinclude/pviewer.h" 446 #include "sys.h" 447 448 static int MatView_SeqDense_ASCII(Mat A,Viewer viewer) 449 { 450 Mat_SeqDense *a = (Mat_SeqDense *) A->data; 451 int ierr, i, j, format; 452 FILE *fd; 453 char *outputname; 454 Scalar *v; 455 456 ierr = ViewerASCIIGetPointer(viewer,&fd); CHKERRQ(ierr); 457 ierr = ViewerFileGetOutputname_Private(viewer,&outputname); CHKERRQ(ierr); 458 ierr = ViewerGetFormat(viewer,&format); 459 if (format == ASCII_FORMAT_INFO || format == ASCII_FORMAT_INFO_DETAILED) { 460 return 0; /* do nothing for now */ 461 } 462 else if (format == ASCII_FORMAT_COMMON) { 463 for ( i=0; i<a->m; i++ ) { 464 v = a->v + i; 465 fprintf(fd,"row %d:",i); 466 for ( j=0; j<a->n; j++ ) { 467 #if defined(PETSC_COMPLEX) 468 if (real(*v) != 0.0 && imag(*v) != 0.0) 469 fprintf(fd," %d %g + %g i",j,real(*v),imag(*v)); 470 else if (real(*v)) fprintf(fd," %d %g ",j,real(*v)); 471 v += a->m; 472 #else 473 if (*v) fprintf(fd," %d %g ",j,*v); 474 v += a->m; 475 #endif 476 } 477 fprintf(fd,"\n"); 478 } 479 } 480 else { 481 #if defined(PETSC_COMPLEX) 482 int allreal = 1; 483 /* determine if matrix has all real values */ 484 v = a->v; 485 for ( i=0; i<a->m*a->n; i++ ) { 486 if (imag(v[i])) { allreal = 0; break ;} 487 } 488 #endif 489 for ( i=0; i<a->m; i++ ) { 490 v = a->v + i; 491 for ( j=0; j<a->n; j++ ) { 492 #if defined(PETSC_COMPLEX) 493 if (allreal) { 494 fprintf(fd,"%6.4e ",real(*v)); v += a->m; 495 } else { 496 fprintf(fd,"%6.4e + %6.4e i ",real(*v),imag(*v)); v += a->m; 497 } 498 #else 499 fprintf(fd,"%6.4e ",*v); v += a->m; 500 #endif 501 } 502 fprintf(fd,"\n"); 503 } 504 } 505 fflush(fd); 506 return 0; 507 } 508 509 static int MatView_SeqDense_Binary(Mat A,Viewer viewer) 510 { 511 Mat_SeqDense *a = (Mat_SeqDense *) A->data; 512 int ict, j, n = a->n, m = a->m, i, fd, *col_lens, ierr, nz = m*n; 513 int format; 514 Scalar *v, *anonz,*vals; 515 516 ierr = ViewerBinaryGetDescriptor(viewer,&fd); CHKERRQ(ierr); 517 518 ierr = ViewerGetFormat(viewer,&format); CHKERRQ(ierr); 519 if (format == BINARY_FORMAT_NATIVE) { 520 /* store the matrix as a dense matrix */ 521 col_lens = (int *) PetscMalloc( 4*sizeof(int) ); CHKPTRQ(col_lens); 522 col_lens[0] = MAT_COOKIE; 523 col_lens[1] = m; 524 col_lens[2] = n; 525 col_lens[3] = MATRIX_BINARY_FORMAT_DENSE; 526 ierr = PetscBinaryWrite(fd,col_lens,4,BINARY_INT,1); CHKERRQ(ierr); 527 PetscFree(col_lens); 528 529 /* write out matrix, by rows */ 530 vals = (Scalar *) PetscMalloc(m*n*sizeof(Scalar)); CHKPTRQ(vals); 531 v = a->v; 532 for ( i=0; i<m; i++ ) { 533 for ( j=0; j<n; j++ ) { 534 vals[i + j*m] = *v++; 535 } 536 } 537 ierr = PetscBinaryWrite(fd,vals,n*m,BINARY_SCALAR,0); CHKERRQ(ierr); 538 PetscFree(vals); 539 } else { 540 col_lens = (int *) PetscMalloc( (4+nz)*sizeof(int) ); CHKPTRQ(col_lens); 541 col_lens[0] = MAT_COOKIE; 542 col_lens[1] = m; 543 col_lens[2] = n; 544 col_lens[3] = nz; 545 546 /* store lengths of each row and write (including header) to file */ 547 for ( i=0; i<m; i++ ) col_lens[4+i] = n; 548 ierr = PetscBinaryWrite(fd,col_lens,4+m,BINARY_INT,1); CHKERRQ(ierr); 549 550 /* Possibly should write in smaller increments, not whole matrix at once? */ 551 /* store column indices (zero start index) */ 552 ict = 0; 553 for ( i=0; i<m; i++ ) { 554 for ( j=0; j<n; j++ ) col_lens[ict++] = j; 555 } 556 ierr = PetscBinaryWrite(fd,col_lens,nz,BINARY_INT,0); CHKERRQ(ierr); 557 PetscFree(col_lens); 558 559 /* store nonzero values */ 560 anonz = (Scalar *) PetscMalloc((nz)*sizeof(Scalar)); CHKPTRQ(anonz); 561 ict = 0; 562 for ( i=0; i<m; i++ ) { 563 v = a->v + i; 564 for ( j=0; j<n; j++ ) { 565 anonz[ict++] = *v; v += a->m; 566 } 567 } 568 ierr = PetscBinaryWrite(fd,anonz,nz,BINARY_SCALAR,0); CHKERRQ(ierr); 569 PetscFree(anonz); 570 } 571 return 0; 572 } 573 574 static int MatView_SeqDense(PetscObject obj,Viewer viewer) 575 { 576 Mat A = (Mat) obj; 577 Mat_SeqDense *a = (Mat_SeqDense*) A->data; 578 ViewerType vtype; 579 int ierr; 580 581 if (!viewer) { 582 viewer = STDOUT_VIEWER_SELF; 583 } 584 ierr = ViewerGetType(viewer,&vtype); CHKERRQ(ierr); 585 586 if (vtype == MATLAB_VIEWER) { 587 return ViewerMatlabPutArray_Private(viewer,a->m,a->n,a->v); 588 } 589 else if (vtype == ASCII_FILE_VIEWER || vtype == ASCII_FILES_VIEWER) { 590 return MatView_SeqDense_ASCII(A,viewer); 591 } 592 else if (vtype == BINARY_FILE_VIEWER) { 593 return MatView_SeqDense_Binary(A,viewer); 594 } 595 return 0; 596 } 597 598 static int MatDestroy_SeqDense(PetscObject obj) 599 { 600 Mat mat = (Mat) obj; 601 Mat_SeqDense *l = (Mat_SeqDense *) mat->data; 602 #if defined(PETSC_LOG) 603 PLogObjectState(obj,"Rows %d Cols %d",l->m,l->n); 604 #endif 605 if (l->pivots) PetscFree(l->pivots); 606 if (!l->user_alloc) PetscFree(l->v); 607 PetscFree(l); 608 PLogObjectDestroy(mat); 609 PetscHeaderDestroy(mat); 610 return 0; 611 } 612 613 static int MatTranspose_SeqDense(Mat A,Mat *matout) 614 { 615 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 616 int k, j, m, n; 617 Scalar *v, tmp; 618 619 v = mat->v; m = mat->m; n = mat->n; 620 if (matout == PETSC_NULL) { /* in place transpose */ 621 if (m != n) SETERRQ(1,"MatTranspose_SeqDense:Supports square matrix only in-place"); 622 for ( j=0; j<m; j++ ) { 623 for ( k=0; k<j; k++ ) { 624 tmp = v[j + k*n]; 625 v[j + k*n] = v[k + j*n]; 626 v[k + j*n] = tmp; 627 } 628 } 629 } 630 else { /* out-of-place transpose */ 631 int ierr; 632 Mat tmat; 633 Mat_SeqDense *tmatd; 634 Scalar *v2; 635 ierr = MatCreateSeqDense(A->comm,mat->n,mat->m,PETSC_NULL,&tmat); CHKERRQ(ierr); 636 tmatd = (Mat_SeqDense *) tmat->data; 637 v = mat->v; v2 = tmatd->v; 638 for ( j=0; j<n; j++ ) { 639 for ( k=0; k<m; k++ ) v2[j + k*n] = v[k + j*m]; 640 } 641 ierr = MatAssemblyBegin(tmat,FINAL_ASSEMBLY); CHKERRQ(ierr); 642 ierr = MatAssemblyEnd(tmat,FINAL_ASSEMBLY); CHKERRQ(ierr); 643 *matout = tmat; 644 } 645 return 0; 646 } 647 648 static int MatEqual_SeqDense(Mat A1,Mat A2, PetscTruth *flg) 649 { 650 Mat_SeqDense *mat1 = (Mat_SeqDense *) A1->data; 651 Mat_SeqDense *mat2 = (Mat_SeqDense *) A2->data; 652 int i; 653 Scalar *v1 = mat1->v, *v2 = mat2->v; 654 655 if (A2->type != MATSEQDENSE) SETERRQ(1,"MatEqual_SeqDense:Both matrices should be of type MATSEQDENSE"); 656 if (mat1->m != mat2->m) {*flg = PETSC_FALSE; return 0;} 657 if (mat1->n != mat2->n) {*flg = PETSC_FALSE; return 0;} 658 for ( i=0; i<mat1->m*mat1->n; i++ ) { 659 if (*v1 != *v2) {*flg = PETSC_FALSE; return 0;} 660 v1++; v2++; 661 } 662 *flg = PETSC_TRUE; 663 return 0; 664 } 665 666 static int MatGetDiagonal_SeqDense(Mat A,Vec v) 667 { 668 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 669 int i, n, len; 670 Scalar *x, zero = 0.0; 671 672 VecSet(&zero,v); 673 VecGetArray(v,&x); VecGetSize(v,&n); 674 len = PetscMin(mat->m,mat->n); 675 if (n != mat->m) SETERRQ(1,"MatGetDiagonal_SeqDense:Nonconforming mat and vec"); 676 for ( i=0; i<len; i++ ) { 677 x[i] = mat->v[i*mat->m + i]; 678 } 679 return 0; 680 } 681 682 static int MatDiagonalScale_SeqDense(Mat A,Vec ll,Vec rr) 683 { 684 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 685 Scalar *l,*r,x,*v; 686 int i,j,m = mat->m, n = mat->n; 687 688 if (ll) { 689 VecGetArray(ll,&l); VecGetSize(ll,&m); 690 if (m != mat->m) SETERRQ(1,"MatDiagonalScale_SeqDense:Left scaling vec wrong size"); 691 for ( i=0; i<m; i++ ) { 692 x = l[i]; 693 v = mat->v + i; 694 for ( j=0; j<n; j++ ) { (*v) *= x; v+= m;} 695 } 696 PLogFlops(n*m); 697 } 698 if (rr) { 699 VecGetArray(rr,&r); VecGetSize(rr,&n); 700 if (n != mat->n) SETERRQ(1,"MatDiagonalScale_SeqDense:Right scaling vec wrong size"); 701 for ( i=0; i<n; i++ ) { 702 x = r[i]; 703 v = mat->v + i*m; 704 for ( j=0; j<m; j++ ) { (*v++) *= x;} 705 } 706 PLogFlops(n*m); 707 } 708 return 0; 709 } 710 711 static int MatNorm_SeqDense(Mat A,NormType type,double *norm) 712 { 713 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 714 Scalar *v = mat->v; 715 double sum = 0.0; 716 int i, j; 717 718 if (type == NORM_FROBENIUS) { 719 for (i=0; i<mat->n*mat->m; i++ ) { 720 #if defined(PETSC_COMPLEX) 721 sum += real(conj(*v)*(*v)); v++; 722 #else 723 sum += (*v)*(*v); v++; 724 #endif 725 } 726 *norm = sqrt(sum); 727 PLogFlops(2*mat->n*mat->m); 728 } 729 else if (type == NORM_1) { 730 *norm = 0.0; 731 for ( j=0; j<mat->n; j++ ) { 732 sum = 0.0; 733 for ( i=0; i<mat->m; i++ ) { 734 sum += PetscAbsScalar(*v); v++; 735 } 736 if (sum > *norm) *norm = sum; 737 } 738 PLogFlops(mat->n*mat->m); 739 } 740 else if (type == NORM_INFINITY) { 741 *norm = 0.0; 742 for ( j=0; j<mat->m; j++ ) { 743 v = mat->v + j; 744 sum = 0.0; 745 for ( i=0; i<mat->n; i++ ) { 746 sum += PetscAbsScalar(*v); v += mat->m; 747 } 748 if (sum > *norm) *norm = sum; 749 } 750 PLogFlops(mat->n*mat->m); 751 } 752 else { 753 SETERRQ(1,"MatNorm_SeqDense:No two norm"); 754 } 755 return 0; 756 } 757 758 static int MatSetOption_SeqDense(Mat A,MatOption op) 759 { 760 Mat_SeqDense *aij = (Mat_SeqDense *) A->data; 761 762 if (op == ROW_ORIENTED) aij->roworiented = 1; 763 else if (op == COLUMN_ORIENTED) aij->roworiented = 0; 764 else if (op == ROWS_SORTED || 765 op == COLUMNS_SORTED || 766 op == SYMMETRIC_MATRIX || 767 op == STRUCTURALLY_SYMMETRIC_MATRIX || 768 op == NO_NEW_NONZERO_LOCATIONS || 769 op == YES_NEW_NONZERO_LOCATIONS || 770 op == NO_NEW_DIAGONALS || 771 op == YES_NEW_DIAGONALS) 772 PLogInfo(A,"Info:MatSetOption_SeqDense:Option ignored\n"); 773 else 774 {SETERRQ(PETSC_ERR_SUP,"MatSetOption_SeqDense:unknown option");} 775 return 0; 776 } 777 778 static int MatZeroEntries_SeqDense(Mat A) 779 { 780 Mat_SeqDense *l = (Mat_SeqDense *) A->data; 781 PetscMemzero(l->v,l->m*l->n*sizeof(Scalar)); 782 return 0; 783 } 784 785 static int MatZeroRows_SeqDense(Mat A,IS is,Scalar *diag) 786 { 787 Mat_SeqDense *l = (Mat_SeqDense *) A->data; 788 int n = l->n, i, j,ierr,N, *rows; 789 Scalar *slot; 790 791 ierr = ISGetSize(is,&N); CHKERRQ(ierr); 792 ierr = ISGetIndices(is,&rows); CHKERRQ(ierr); 793 for ( i=0; i<N; i++ ) { 794 slot = l->v + rows[i]; 795 for ( j=0; j<n; j++ ) { *slot = 0.0; slot += n;} 796 } 797 if (diag) { 798 for ( i=0; i<N; i++ ) { 799 slot = l->v + (n+1)*rows[i]; 800 *slot = *diag; 801 } 802 } 803 ISRestoreIndices(is,&rows); 804 return 0; 805 } 806 807 static int MatGetSize_SeqDense(Mat A,int *m,int *n) 808 { 809 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 810 *m = mat->m; *n = mat->n; 811 return 0; 812 } 813 814 static int MatGetOwnershipRange_SeqDense(Mat A,int *m,int *n) 815 { 816 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 817 *m = 0; *n = mat->m; 818 return 0; 819 } 820 821 static int MatGetArray_SeqDense(Mat A,Scalar **array) 822 { 823 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 824 *array = mat->v; 825 return 0; 826 } 827 828 static int MatRestoreArray_SeqDense(Mat A,Scalar **array) 829 { 830 return 0; 831 } 832 833 static int MatGetSubMatrixInPlace_SeqDense(Mat A,IS isrow,IS iscol) 834 { 835 SETERRQ(1,"MatGetSubMatrixInPlace_SeqDense:not done"); 836 } 837 838 static int MatGetSubMatrix_SeqDense(Mat A,IS isrow,IS iscol,MatGetSubMatrixCall scall, 839 Mat *submat) 840 { 841 Mat_SeqDense *mat = (Mat_SeqDense *) A->data; 842 int nznew, *smap, i, j, ierr, oldcols = mat->n; 843 int *irow, *icol, nrows, ncols, *cwork; 844 Scalar *vwork, *val; 845 Mat newmat; 846 847 ierr = ISGetIndices(isrow,&irow); CHKERRQ(ierr); 848 ierr = ISGetIndices(iscol,&icol); CHKERRQ(ierr); 849 ierr = ISGetSize(isrow,&nrows); CHKERRQ(ierr); 850 ierr = ISGetSize(iscol,&ncols); CHKERRQ(ierr); 851 852 smap = (int *) PetscMalloc(oldcols*sizeof(int)); CHKPTRQ(smap); 853 cwork = (int *) PetscMalloc(ncols*sizeof(int)); CHKPTRQ(cwork); 854 vwork = (Scalar *) PetscMalloc(ncols*sizeof(Scalar)); CHKPTRQ(vwork); 855 PetscMemzero((char*)smap,oldcols*sizeof(int)); 856 for ( i=0; i<ncols; i++ ) smap[icol[i]] = i+1; 857 858 /* Create and fill new matrix */ 859 ierr = MatCreateSeqDense(A->comm,nrows,ncols,PETSC_NULL,&newmat); CHKERRQ(ierr); 860 for (i=0; i<nrows; i++) { 861 nznew = 0; 862 val = mat->v + irow[i]; 863 for (j=0; j<oldcols; j++) { 864 if (smap[j]) { 865 cwork[nznew] = smap[j] - 1; 866 vwork[nznew++] = val[j * mat->m]; 867 } 868 } 869 ierr = MatSetValues(newmat,1,&i,nznew,cwork,vwork,INSERT_VALUES); 870 CHKERRQ(ierr); 871 } 872 ierr = MatAssemblyBegin(newmat,FINAL_ASSEMBLY); CHKERRQ(ierr); 873 ierr = MatAssemblyEnd(newmat,FINAL_ASSEMBLY); CHKERRQ(ierr); 874 875 /* Free work space */ 876 PetscFree(smap); PetscFree(cwork); PetscFree(vwork); 877 ierr = ISRestoreIndices(isrow,&irow); CHKERRQ(ierr); 878 ierr = ISRestoreIndices(iscol,&icol); CHKERRQ(ierr); 879 *submat = newmat; 880 return 0; 881 } 882 883 static int MatCopy_SeqDense(Mat A, Mat B) 884 { 885 Mat_SeqDense *a = (Mat_SeqDense *) A->data, *b = (Mat_SeqDense *)B->data; 886 if (B->type != MATSEQDENSE) return MatCopy_Basic(A,B); 887 if (a->m != b->m || a->n != b->n) SETERRQ(1,"MatCopy_SeqDense:size(B) != size(A)"); 888 PetscMemcpy(b->v,a->v,a->m*a->n*sizeof(Scalar)); 889 return 0; 890 } 891 892 /* -------------------------------------------------------------------*/ 893 static struct _MatOps MatOps = {MatSetValues_SeqDense, 894 MatGetRow_SeqDense, MatRestoreRow_SeqDense, 895 MatMult_SeqDense, MatMultAdd_SeqDense, 896 MatMultTrans_SeqDense, MatMultTransAdd_SeqDense, 897 MatSolve_SeqDense,MatSolveAdd_SeqDense, 898 MatSolveTrans_SeqDense,MatSolveTransAdd_SeqDense, 899 MatLUFactor_SeqDense,MatCholeskyFactor_SeqDense, 900 MatRelax_SeqDense, 901 MatTranspose_SeqDense, 902 MatGetInfo_SeqDense,MatEqual_SeqDense, 903 MatGetDiagonal_SeqDense,MatDiagonalScale_SeqDense,MatNorm_SeqDense, 904 0,0, 905 0, MatSetOption_SeqDense,MatZeroEntries_SeqDense,MatZeroRows_SeqDense,0, 906 MatLUFactorSymbolic_SeqDense,MatLUFactorNumeric_SeqDense, 907 MatCholeskyFactorSymbolic_SeqDense,MatCholeskyFactorNumeric_SeqDense, 908 MatGetSize_SeqDense,MatGetSize_SeqDense,MatGetOwnershipRange_SeqDense, 909 0,0, MatGetArray_SeqDense, MatRestoreArray_SeqDense,0, 910 MatGetSubMatrix_SeqDense,MatGetSubMatrixInPlace_SeqDense, 911 MatConvertSameType_SeqDense,0,0,0,0, 912 MatAXPY_SeqDense,0,0, 913 MatGetValues_SeqDense, 914 MatCopy_SeqDense,0,MatScale_SeqDense}; 915 916 917 /*@C 918 MatCreateSeqDense - Creates a sequential dense matrix that 919 is stored in column major order (the usual Fortran 77 manner). Many 920 of the matrix operations use the BLAS and LAPACK routines. 921 922 Input Parameters: 923 . comm - MPI communicator, set to MPI_COMM_SELF 924 . m - number of rows 925 . n - number of columns 926 . data - optional location of matrix data. Set data=PETSC_NULL for PETSc 927 to control all matrix memory allocation. 928 929 Output Parameter: 930 . A - the matrix 931 932 Notes: 933 The data input variable is intended primarily for Fortran programmers 934 who wish to allocate their own matrix memory space. Most users should 935 set data=PETSC_NULL. 936 937 .keywords: dense, matrix, LAPACK, BLAS 938 939 .seealso: MatCreate(), MatSetValues() 940 @*/ 941 int MatCreateSeqDense(MPI_Comm comm,int m,int n,Scalar *data,Mat *A) 942 { 943 Mat B; 944 Mat_SeqDense *b; 945 int ierr,flg; 946 947 *A = 0; 948 PetscHeaderCreate(B,_Mat,MAT_COOKIE,MATSEQDENSE,comm); 949 PLogObjectCreate(B); 950 b = (Mat_SeqDense *) PetscMalloc(sizeof(Mat_SeqDense)); CHKPTRQ(b); 951 PetscMemzero(b,sizeof(Mat_SeqDense)); 952 PetscMemcpy(&B->ops,&MatOps,sizeof(struct _MatOps)); 953 B->destroy = MatDestroy_SeqDense; 954 B->view = MatView_SeqDense; 955 B->factor = 0; 956 PLogObjectMemory(B,sizeof(struct _Mat)); 957 B->data = (void *) b; 958 959 b->m = m; B->m = m; B->M = m; 960 b->n = n; B->n = n; B->N = n; 961 b->pivots = 0; 962 b->roworiented = 1; 963 if (data == PETSC_NULL) { 964 b->v = (Scalar*) PetscMalloc((m*n+1)*sizeof(Scalar)); CHKPTRQ(b->v); 965 PetscMemzero(b->v,m*n*sizeof(Scalar)); 966 b->user_alloc = 0; 967 PLogObjectMemory(B,n*m*sizeof(Scalar)); 968 } 969 else { /* user-allocated storage */ 970 b->v = data; 971 b->user_alloc = 1; 972 } 973 ierr = OptionsHasName(PETSC_NULL,"-help",&flg); CHKERRQ(ierr); 974 if (flg) {ierr = MatPrintHelp(B); CHKERRQ(ierr);} 975 *A = B; 976 return 0; 977 } 978 979 int MatCreate_SeqDense(Mat A,Mat *newmat) 980 { 981 Mat_SeqDense *m = (Mat_SeqDense *) A->data; 982 return MatCreateSeqDense(A->comm,m->m,m->n,PETSC_NULL,newmat); 983 } 984