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