1 #ifndef lint 2 static char vcid[] = "$Id: aij.c,v 1.222 1997/06/18 01:44:15 curfman Exp curfman $"; 3 #endif 4 5 /* 6 Defines the basic matrix operations for the AIJ (compressed row) 7 matrix storage format. 8 */ 9 10 #include "pinclude/pviewer.h" 11 #include "sys.h" 12 #include "src/mat/impls/aij/seq/aij.h" 13 #include "src/vec/vecimpl.h" 14 #include "src/inline/spops.h" 15 #include "petsc.h" 16 #include "src/inline/bitarray.h" 17 18 /* 19 Basic AIJ format ILU based on drop tolerance 20 */ 21 #undef __FUNC__ 22 #define __FUNC__ "MatILUDTFactor_SeqAIJ" 23 int MatILUDTFactor_SeqAIJ(Mat A,double dt,int maxnz,IS row,IS col,Mat *fact) 24 { 25 /* Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; */ 26 int ierr = 1; 27 28 SETERRQ(ierr,0,"Not implemented"); 29 } 30 31 extern int MatToSymmetricIJ_SeqAIJ(int,int*,int*,int,int,int**,int**); 32 33 #undef __FUNC__ 34 #define __FUNC__ "MatGetRowIJ_SeqAIJ" 35 int MatGetRowIJ_SeqAIJ(Mat A,int oshift,PetscTruth symmetric,int *m,int **ia,int **ja, 36 PetscTruth *done) 37 { 38 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 39 int ierr,i,ishift; 40 41 *m = A->m; 42 if (!ia) return 0; 43 ishift = a->indexshift; 44 if (symmetric) { 45 ierr = MatToSymmetricIJ_SeqAIJ(a->m,a->i,a->j,ishift,oshift,ia,ja); CHKERRQ(ierr); 46 } else if (oshift == 0 && ishift == -1) { 47 int nz = a->i[a->m]; 48 /* malloc space and subtract 1 from i and j indices */ 49 *ia = (int *) PetscMalloc( (a->m+1)*sizeof(int) ); CHKPTRQ(*ia); 50 *ja = (int *) PetscMalloc( (nz+1)*sizeof(int) ); CHKPTRQ(*ja); 51 for ( i=0; i<nz; i++ ) (*ja)[i] = a->j[i] - 1; 52 for ( i=0; i<a->m+1; i++ ) (*ia)[i] = a->i[i] - 1; 53 } else if (oshift == 1 && ishift == 0) { 54 int nz = a->i[a->m] + 1; 55 /* malloc space and add 1 to i and j indices */ 56 *ia = (int *) PetscMalloc( (a->m+1)*sizeof(int) ); CHKPTRQ(*ia); 57 *ja = (int *) PetscMalloc( (nz+1)*sizeof(int) ); CHKPTRQ(*ja); 58 for ( i=0; i<nz; i++ ) (*ja)[i] = a->j[i] + 1; 59 for ( i=0; i<a->m+1; i++ ) (*ia)[i] = a->i[i] + 1; 60 } else { 61 *ia = a->i; *ja = a->j; 62 } 63 64 return 0; 65 } 66 67 #undef __FUNC__ 68 #define __FUNC__ "MatRestoreRowIJ_SeqAIJ" 69 int MatRestoreRowIJ_SeqAIJ(Mat A,int oshift,PetscTruth symmetric,int *n,int **ia,int **ja, 70 PetscTruth *done) 71 { 72 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 73 int ishift = a->indexshift; 74 75 if (!ia) return 0; 76 if (symmetric || (oshift == 0 && ishift == -1) || (oshift == 1 && ishift == 0)) { 77 PetscFree(*ia); 78 PetscFree(*ja); 79 } 80 return 0; 81 } 82 83 #undef __FUNC__ 84 #define __FUNC__ "MatGetColumnIJ_SeqAIJ" 85 int MatGetColumnIJ_SeqAIJ(Mat A,int oshift,PetscTruth symmetric,int *nn,int **ia,int **ja, 86 PetscTruth *done) 87 { 88 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 89 int ierr,i,ishift = a->indexshift,*collengths,*cia,*cja,n = A->n,m = A->m; 90 int nz = a->i[m]+ishift,row,*jj,mr,col; 91 92 *nn = A->n; 93 if (!ia) return 0; 94 if (symmetric) { 95 ierr = MatToSymmetricIJ_SeqAIJ(a->m,a->i,a->j,ishift,oshift,ia,ja); CHKERRQ(ierr); 96 } else { 97 collengths = (int *) PetscMalloc( (n+1)*sizeof(int) ); CHKPTRQ(collengths); 98 PetscMemzero(collengths,n*sizeof(int)); 99 cia = (int *) PetscMalloc( (n+1)*sizeof(int) ); CHKPTRQ(cia); 100 cja = (int *) PetscMalloc( (nz+1)*sizeof(int) ); CHKPTRQ(cja); 101 jj = a->j; 102 for ( i=0; i<nz; i++ ) { 103 collengths[jj[i] + ishift]++; 104 } 105 cia[0] = oshift; 106 for ( i=0; i<n; i++) { 107 cia[i+1] = cia[i] + collengths[i]; 108 } 109 PetscMemzero(collengths,n*sizeof(int)); 110 jj = a->j; 111 for ( row=0; row<m; row++ ) { 112 mr = a->i[row+1] - a->i[row]; 113 for ( i=0; i<mr; i++ ) { 114 col = *jj++ + ishift; 115 cja[cia[col] + collengths[col]++ - oshift] = row + oshift; 116 } 117 } 118 PetscFree(collengths); 119 *ia = cia; *ja = cja; 120 } 121 122 return 0; 123 } 124 125 #undef __FUNC__ 126 #define __FUNC__ "MatRestoreColumnIJ_SeqAIJ" 127 int MatRestoreColumnIJ_SeqAIJ(Mat A,int oshift,PetscTruth symmetric,int *n,int **ia, 128 int **ja,PetscTruth *done) 129 { 130 if (!ia) return 0; 131 132 PetscFree(*ia); 133 PetscFree(*ja); 134 135 return 0; 136 } 137 138 #define CHUNKSIZE 15 139 140 #undef __FUNC__ 141 #define __FUNC__ "MatSetValues_SeqAIJ" 142 int MatSetValues_SeqAIJ(Mat A,int m,int *im,int n,int *in,Scalar *v,InsertMode is) 143 { 144 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 145 int *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax, N, sorted = a->sorted; 146 int *imax = a->imax, *ai = a->i, *ailen = a->ilen,roworiented = a->roworiented; 147 int *aj = a->j, nonew = a->nonew,shift = a->indexshift; 148 Scalar *ap,value, *aa = a->a; 149 150 for ( k=0; k<m; k++ ) { /* loop over added rows */ 151 row = im[k]; 152 #if defined(PETSC_BOPT_g) 153 if (row < 0) SETERRQ(1,0,"Negative row"); 154 if (row >= a->m) SETERRQ(1,0,"Row too large"); 155 #endif 156 rp = aj + ai[row] + shift; ap = aa + ai[row] + shift; 157 rmax = imax[row]; nrow = ailen[row]; 158 low = 0; 159 for ( l=0; l<n; l++ ) { /* loop over added columns */ 160 #if defined(PETSC_BOPT_g) 161 if (in[l] < 0) SETERRQ(1,0,"Negative column"); 162 if (in[l] >= a->n) SETERRQ(1,0,"Column too large"); 163 #endif 164 col = in[l] - shift; 165 if (roworiented) { 166 value = *v++; 167 } 168 else { 169 value = v[k + l*m]; 170 } 171 if (!sorted) low = 0; high = nrow; 172 while (high-low > 5) { 173 t = (low+high)/2; 174 if (rp[t] > col) high = t; 175 else low = t; 176 } 177 for ( i=low; i<high; i++ ) { 178 if (rp[i] > col) break; 179 if (rp[i] == col) { 180 if (is == ADD_VALUES) ap[i] += value; 181 else ap[i] = value; 182 goto noinsert; 183 } 184 } 185 if (nonew == 1) goto noinsert; 186 else if (nonew == -1) SETERRQ(1,0,"Inserting a new nonzero in the matrix"); 187 if (nrow >= rmax) { 188 /* there is no extra room in row, therefore enlarge */ 189 int new_nz = ai[a->m] + CHUNKSIZE,len,*new_i,*new_j; 190 Scalar *new_a; 191 192 if (nonew == -2) SETERRQ(1,0,"Inserting a new nonzero in the matrix"); 193 194 /* malloc new storage space */ 195 len = new_nz*(sizeof(int)+sizeof(Scalar))+(a->m+1)*sizeof(int); 196 new_a = (Scalar *) PetscMalloc( len ); CHKPTRQ(new_a); 197 new_j = (int *) (new_a + new_nz); 198 new_i = new_j + new_nz; 199 200 /* copy over old data into new slots */ 201 for ( ii=0; ii<row+1; ii++ ) {new_i[ii] = ai[ii];} 202 for ( ii=row+1; ii<a->m+1; ii++ ) {new_i[ii] = ai[ii]+CHUNKSIZE;} 203 PetscMemcpy(new_j,aj,(ai[row]+nrow+shift)*sizeof(int)); 204 len = (new_nz - CHUNKSIZE - ai[row] - nrow - shift); 205 PetscMemcpy(new_j+ai[row]+shift+nrow+CHUNKSIZE,aj+ai[row]+shift+nrow, 206 len*sizeof(int)); 207 PetscMemcpy(new_a,aa,(ai[row]+nrow+shift)*sizeof(Scalar)); 208 PetscMemcpy(new_a+ai[row]+shift+nrow+CHUNKSIZE,aa+ai[row]+shift+nrow, 209 len*sizeof(Scalar)); 210 /* free up old matrix storage */ 211 PetscFree(a->a); 212 if (!a->singlemalloc) {PetscFree(a->i);PetscFree(a->j);} 213 aa = a->a = new_a; ai = a->i = new_i; aj = a->j = new_j; 214 a->singlemalloc = 1; 215 216 rp = aj + ai[row] + shift; ap = aa + ai[row] + shift; 217 rmax = imax[row] = imax[row] + CHUNKSIZE; 218 PLogObjectMemory(A,CHUNKSIZE*(sizeof(int) + sizeof(Scalar))); 219 a->maxnz += CHUNKSIZE; 220 a->reallocs++; 221 } 222 N = nrow++ - 1; a->nz++; 223 /* shift up all the later entries in this row */ 224 for ( ii=N; ii>=i; ii-- ) { 225 rp[ii+1] = rp[ii]; 226 ap[ii+1] = ap[ii]; 227 } 228 rp[i] = col; 229 ap[i] = value; 230 noinsert:; 231 low = i + 1; 232 } 233 ailen[row] = nrow; 234 } 235 return 0; 236 } 237 238 #undef __FUNC__ 239 #define __FUNC__ "MatGetValues_SeqAIJ" 240 int MatGetValues_SeqAIJ(Mat A,int m,int *im,int n,int *in,Scalar *v) 241 { 242 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 243 int *rp, k, low, high, t, row, nrow, i, col, l, *aj = a->j; 244 int *ai = a->i, *ailen = a->ilen, shift = a->indexshift; 245 Scalar *ap, *aa = a->a, zero = 0.0; 246 247 for ( k=0; k<m; k++ ) { /* loop over rows */ 248 row = im[k]; 249 if (row < 0) SETERRQ(1,0,"Negative row"); 250 if (row >= a->m) SETERRQ(1,0,"Row too large"); 251 rp = aj + ai[row] + shift; ap = aa + ai[row] + shift; 252 nrow = ailen[row]; 253 for ( l=0; l<n; l++ ) { /* loop over columns */ 254 if (in[l] < 0) SETERRQ(1,0,"Negative column"); 255 if (in[l] >= a->n) SETERRQ(1,0,"Column too large"); 256 col = in[l] - shift; 257 high = nrow; low = 0; /* assume unsorted */ 258 while (high-low > 5) { 259 t = (low+high)/2; 260 if (rp[t] > col) high = t; 261 else low = t; 262 } 263 for ( i=low; i<high; i++ ) { 264 if (rp[i] > col) break; 265 if (rp[i] == col) { 266 *v++ = ap[i]; 267 goto finished; 268 } 269 } 270 *v++ = zero; 271 finished:; 272 } 273 } 274 return 0; 275 } 276 277 278 #undef __FUNC__ 279 #define __FUNC__ "MatView_SeqAIJ_Binary" 280 extern int MatView_SeqAIJ_Binary(Mat A,Viewer viewer) 281 { 282 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 283 int i, fd, *col_lens, ierr; 284 285 ierr = ViewerBinaryGetDescriptor(viewer,&fd); CHKERRQ(ierr); 286 col_lens = (int *) PetscMalloc( (4+a->m)*sizeof(int) ); CHKPTRQ(col_lens); 287 col_lens[0] = MAT_COOKIE; 288 col_lens[1] = a->m; 289 col_lens[2] = a->n; 290 col_lens[3] = a->nz; 291 292 /* store lengths of each row and write (including header) to file */ 293 for ( i=0; i<a->m; i++ ) { 294 col_lens[4+i] = a->i[i+1] - a->i[i]; 295 } 296 ierr = PetscBinaryWrite(fd,col_lens,4+a->m,BINARY_INT,1); CHKERRQ(ierr); 297 PetscFree(col_lens); 298 299 /* store column indices (zero start index) */ 300 if (a->indexshift) { 301 for ( i=0; i<a->nz; i++ ) a->j[i]--; 302 } 303 ierr = PetscBinaryWrite(fd,a->j,a->nz,BINARY_INT,0); CHKERRQ(ierr); 304 if (a->indexshift) { 305 for ( i=0; i<a->nz; i++ ) a->j[i]++; 306 } 307 308 /* store nonzero values */ 309 ierr = PetscBinaryWrite(fd,a->a,a->nz,BINARY_SCALAR,0); CHKERRQ(ierr); 310 return 0; 311 } 312 313 #undef __FUNC__ 314 #define __FUNC__ "MatView_SeqAIJ_ASCII" 315 extern int MatView_SeqAIJ_ASCII(Mat A,Viewer viewer) 316 { 317 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 318 int ierr, i,j, m = a->m, shift = a->indexshift, format, flg1,flg2; 319 FILE *fd; 320 char *outputname; 321 322 ierr = ViewerASCIIGetPointer(viewer,&fd); CHKERRQ(ierr); 323 ierr = ViewerFileGetOutputname_Private(viewer,&outputname); CHKERRQ(ierr); 324 ierr = ViewerGetFormat(viewer,&format); 325 if (format == VIEWER_FORMAT_ASCII_INFO) { 326 return 0; 327 } 328 else if (format == VIEWER_FORMAT_ASCII_INFO_LONG) { 329 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_no_inode",&flg1); CHKERRQ(ierr); 330 ierr = OptionsHasName(PETSC_NULL,"-mat_no_unroll",&flg2); CHKERRQ(ierr); 331 if (flg1 || flg2) fprintf(fd," not using I-node routines\n"); 332 else fprintf(fd," using I-node routines: found %d nodes, limit used is %d\n", 333 a->inode.node_count,a->inode.limit); 334 } 335 else if (format == VIEWER_FORMAT_ASCII_MATLAB) { 336 fprintf(fd,"%% Size = %d %d \n",m,a->n); 337 fprintf(fd,"%% Nonzeros = %d \n",a->nz); 338 fprintf(fd,"zzz = zeros(%d,3);\n",a->nz); 339 fprintf(fd,"zzz = [\n"); 340 341 for (i=0; i<m; i++) { 342 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 343 #if defined(PETSC_COMPLEX) 344 fprintf(fd,"%d %d %18.16e + %18.16e i \n",i+1,a->j[j]+!shift,real(a->a[j]), 345 imag(a->a[j])); 346 #else 347 fprintf(fd,"%d %d %18.16e\n", i+1, a->j[j]+!shift, a->a[j]); 348 #endif 349 } 350 } 351 fprintf(fd,"];\n %s = spconvert(zzz);\n",outputname); 352 } 353 else if (format == VIEWER_FORMAT_ASCII_COMMON) { 354 for ( i=0; i<m; i++ ) { 355 fprintf(fd,"row %d:",i); 356 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 357 #if defined(PETSC_COMPLEX) 358 if (imag(a->a[j]) > 0.0 && real(a->a[j]) != 0.0) 359 fprintf(fd," %d %g + %g i",a->j[j]+shift,real(a->a[j]),imag(a->a[j])); 360 else if (imag(a->a[j]) < 0.0 && real(a->a[j]) != 0.0) 361 fprintf(fd," %d %g - %g i",a->j[j]+shift,real(a->a[j]),-imag(a->a[j])); 362 else if (real(a->a[j]) != 0.0) 363 fprintf(fd," %d %g ",a->j[j]+shift,real(a->a[j])); 364 #else 365 if (a->a[j] != 0.0) fprintf(fd," %d %g ",a->j[j]+shift,a->a[j]); 366 #endif 367 } 368 fprintf(fd,"\n"); 369 } 370 } 371 else if (format == VIEWER_FORMAT_ASCII_SYMMODU) { 372 int nzd=0, fshift=1, *sptr; 373 sptr = (int *) PetscMalloc( (m)*sizeof(int) ); CHKPTRQ(sptr); 374 for ( i=0; i<m; i++ ) { 375 sptr[i] = nzd+1; 376 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 377 if (a->j[j] >= i) { 378 #if defined(PETSC_COMPLEX) 379 if (imag(a->a[j]) != 0.0 || real(a->a[j]) != 0.0) nzd++; 380 #else 381 if (a->a[j] != 0.0) nzd++; 382 #endif 383 } 384 } 385 } 386 fprintf(fd," %d %d\n\n",m,nzd); 387 for ( i=0; i<m; i+=6 ) { 388 fprintf(fd," %d %d %d %d %d %d\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3],sptr[i+4],sptr[i+5]); 389 } 390 fprintf(fd,"\n"); 391 PetscFree(sptr); 392 for ( i=0; i<m; i++ ) { 393 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 394 if (a->j[j] >= i) fprintf(fd," %d ",a->j[j]+fshift); 395 } 396 fprintf(fd,"\n"); 397 } 398 fprintf(fd,"\n"); 399 for ( i=0; i<m; i++ ) { 400 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 401 if (a->j[j] >= i) { 402 #if defined(PETSC_COMPLEX) 403 if (imag(a->a[j]) != 0.0 || real(a->a[j]) != 0.0) 404 fprintf(fd," %18.16e %18.16e ",real(a->a[j]),imag(a->a[j])); 405 #else 406 if (a->a[j] != 0.0) fprintf(fd," %18.16e ",a->a[j]); 407 #endif 408 } 409 } 410 fprintf(fd,"\n"); 411 } 412 } 413 else { 414 for ( i=0; i<m; i++ ) { 415 fprintf(fd,"row %d:",i); 416 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 417 #if defined(PETSC_COMPLEX) 418 if (imag(a->a[j]) > 0.0) { 419 fprintf(fd," %d %g + %g i",a->j[j]+shift,real(a->a[j]),imag(a->a[j])); 420 } else if (imag(a->a[j]) < 0.0) { 421 fprintf(fd," %d %g - %g i",a->j[j]+shift,real(a->a[j]),-imag(a->a[j])); 422 } 423 else { 424 fprintf(fd," %d %g ",a->j[j]+shift,real(a->a[j])); 425 } 426 #else 427 fprintf(fd," %d %g ",a->j[j]+shift,a->a[j]); 428 #endif 429 } 430 fprintf(fd,"\n"); 431 } 432 } 433 fflush(fd); 434 return 0; 435 } 436 437 #undef __FUNC__ 438 #define __FUNC__ "MatView_SeqAIJ_Draw" 439 extern int MatView_SeqAIJ_Draw(Mat A,Viewer viewer) 440 { 441 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 442 int ierr, i,j, m = a->m, shift = a->indexshift,pause,color; 443 int format; 444 double xl,yl,xr,yr,w,h,xc,yc,scale = 1.0,x_l,x_r,y_l,y_r,maxv = 0.0; 445 Draw draw; 446 DrawButton button; 447 PetscTruth isnull; 448 449 ierr = ViewerDrawGetDraw(viewer,&draw); CHKERRQ(ierr); 450 ierr = DrawClear(draw); CHKERRQ(ierr); 451 ierr = ViewerGetFormat(viewer,&format); CHKERRQ(ierr); 452 ierr = DrawIsNull(draw,&isnull); CHKERRQ(ierr); if (isnull) return 0; 453 454 xr = a->n; yr = a->m; h = yr/10.0; w = xr/10.0; 455 xr += w; yr += h; xl = -w; yl = -h; 456 ierr = DrawSetCoordinates(draw,xl,yl,xr,yr); CHKERRQ(ierr); 457 /* loop over matrix elements drawing boxes */ 458 459 if (format != VIEWER_FORMAT_DRAW_CONTOUR) { 460 /* Blue for negative, Cyan for zero and Red for positive */ 461 color = DRAW_BLUE; 462 for ( i=0; i<m; i++ ) { 463 y_l = m - i - 1.0; y_r = y_l + 1.0; 464 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 465 x_l = a->j[j] + shift; x_r = x_l + 1.0; 466 #if defined(PETSC_COMPLEX) 467 if (real(a->a[j]) >= 0.) continue; 468 #else 469 if (a->a[j] >= 0.) continue; 470 #endif 471 DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color); 472 } 473 } 474 color = DRAW_CYAN; 475 for ( i=0; i<m; i++ ) { 476 y_l = m - i - 1.0; y_r = y_l + 1.0; 477 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 478 x_l = a->j[j] + shift; x_r = x_l + 1.0; 479 if (a->a[j] != 0.) continue; 480 DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color); 481 } 482 } 483 color = DRAW_RED; 484 for ( i=0; i<m; i++ ) { 485 y_l = m - i - 1.0; y_r = y_l + 1.0; 486 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 487 x_l = a->j[j] + shift; x_r = x_l + 1.0; 488 #if defined(PETSC_COMPLEX) 489 if (real(a->a[j]) <= 0.) continue; 490 #else 491 if (a->a[j] <= 0.) continue; 492 #endif 493 DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color); 494 } 495 } 496 } else { 497 /* use contour shading to indicate magnitude of values */ 498 /* first determine max of all nonzero values */ 499 int nz = a->nz,count; 500 Draw popup; 501 502 for ( i=0; i<nz; i++ ) { 503 if (PetscAbsScalar(a->a[i]) > maxv) maxv = PetscAbsScalar(a->a[i]); 504 } 505 ierr = DrawCreatePopUp(draw,&popup); CHKERRQ(ierr); 506 ierr = DrawScalePopup(popup,0.0,maxv); CHKERRQ(ierr); 507 count = 0; 508 for ( i=0; i<m; i++ ) { 509 y_l = m - i - 1.0; y_r = y_l + 1.0; 510 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 511 x_l = a->j[j] + shift; x_r = x_l + 1.0; 512 color = 32 + (int) ((200.0 - 32.0)*PetscAbsScalar(a->a[count])/maxv); 513 DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color); 514 count++; 515 } 516 } 517 } 518 DrawFlush(draw); 519 DrawGetPause(draw,&pause); 520 if (pause >= 0) { PetscSleep(pause); return 0;} 521 522 /* allow the matrix to zoom or shrink */ 523 ierr = DrawCheckResizedWindow(draw); 524 ierr = DrawGetMouseButton(draw,&button,&xc,&yc,0,0); 525 while (button != BUTTON_RIGHT) { 526 DrawClear(draw); 527 if (button == BUTTON_LEFT) scale = .5; 528 else if (button == BUTTON_CENTER) scale = 2.; 529 xl = scale*(xl + w - xc) + xc - w*scale; 530 xr = scale*(xr - w - xc) + xc + w*scale; 531 yl = scale*(yl + h - yc) + yc - h*scale; 532 yr = scale*(yr - h - yc) + yc + h*scale; 533 w *= scale; h *= scale; 534 ierr = DrawSetCoordinates(draw,xl,yl,xr,yr); CHKERRQ(ierr); 535 if (format != VIEWER_FORMAT_DRAW_CONTOUR) { 536 /* Blue for negative, Cyan for zero and Red for positive */ 537 color = DRAW_BLUE; 538 for ( i=0; i<m; i++ ) { 539 y_l = m - i - 1.0; y_r = y_l + 1.0; 540 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 541 x_l = a->j[j] + shift; x_r = x_l + 1.0; 542 #if defined(PETSC_COMPLEX) 543 if (real(a->a[j]) >= 0.) continue; 544 #else 545 if (a->a[j] >= 0.) continue; 546 #endif 547 DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color); 548 } 549 } 550 color = DRAW_CYAN; 551 for ( i=0; i<m; i++ ) { 552 y_l = m - i - 1.0; y_r = y_l + 1.0; 553 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 554 x_l = a->j[j] + shift; x_r = x_l + 1.0; 555 if (a->a[j] != 0.) continue; 556 DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color); 557 } 558 } 559 color = DRAW_RED; 560 for ( i=0; i<m; i++ ) { 561 y_l = m - i - 1.0; y_r = y_l + 1.0; 562 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 563 x_l = a->j[j] + shift; x_r = x_l + 1.0; 564 #if defined(PETSC_COMPLEX) 565 if (real(a->a[j]) <= 0.) continue; 566 #else 567 if (a->a[j] <= 0.) continue; 568 #endif 569 DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color); 570 } 571 } 572 } else { 573 /* use contour shading to indicate magnitude of values */ 574 int count = 0; 575 for ( i=0; i<m; i++ ) { 576 y_l = m - i - 1.0; y_r = y_l + 1.0; 577 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 578 x_l = a->j[j] + shift; x_r = x_l + 1.0; 579 color = 32 + (int) ((200.0 - 32.0)*PetscAbsScalar(a->a[count])/maxv); 580 DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color); 581 count++; 582 } 583 } 584 } 585 586 ierr = DrawCheckResizedWindow(draw); 587 ierr = DrawGetMouseButton(draw,&button,&xc,&yc,0,0); 588 } 589 return 0; 590 } 591 592 #undef __FUNC__ 593 #define __FUNC__ "MatView_SeqAIJ" /* ADIC Ignore */ 594 int MatView_SeqAIJ(PetscObject obj,Viewer viewer) 595 { 596 Mat A = (Mat) obj; 597 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 598 ViewerType vtype; 599 int ierr; 600 601 ierr = ViewerGetType(viewer,&vtype); CHKERRQ(ierr); 602 if (vtype == MATLAB_VIEWER) { 603 return ViewerMatlabPutSparse_Private(viewer,a->m,a->n,a->nz,a->a,a->i,a->j); 604 } 605 else if (vtype == ASCII_FILE_VIEWER || vtype == ASCII_FILES_VIEWER){ 606 return MatView_SeqAIJ_ASCII(A,viewer); 607 } 608 else if (vtype == BINARY_FILE_VIEWER) { 609 return MatView_SeqAIJ_Binary(A,viewer); 610 } 611 else if (vtype == DRAW_VIEWER) { 612 return MatView_SeqAIJ_Draw(A,viewer); 613 } 614 return 0; 615 } 616 617 extern int Mat_AIJ_CheckInode(Mat); 618 #undef __FUNC__ 619 #define __FUNC__ "MatAssemblyEnd_SeqAIJ" 620 int MatAssemblyEnd_SeqAIJ(Mat A,MatAssemblyType mode) 621 { 622 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 623 int fshift = 0,i,j,*ai = a->i, *aj = a->j, *imax = a->imax,ierr; 624 int m = a->m, *ip, N, *ailen = a->ilen,shift = a->indexshift,rmax; 625 Scalar *aa = a->a, *ap; 626 627 if (mode == MAT_FLUSH_ASSEMBLY) return 0; 628 629 rmax = ailen[0]; /* determine row with most nonzeros */ 630 for ( i=1; i<m; i++ ) { 631 /* move each row back by the amount of empty slots (fshift) before it*/ 632 fshift += imax[i-1] - ailen[i-1]; 633 rmax = PetscMax(rmax,ailen[i]); 634 if (fshift) { 635 ip = aj + ai[i] + shift; ap = aa + ai[i] + shift; 636 N = ailen[i]; 637 for ( j=0; j<N; j++ ) { 638 ip[j-fshift] = ip[j]; 639 ap[j-fshift] = ap[j]; 640 } 641 } 642 ai[i] = ai[i-1] + ailen[i-1]; 643 } 644 if (m) { 645 fshift += imax[m-1] - ailen[m-1]; 646 ai[m] = ai[m-1] + ailen[m-1]; 647 } 648 /* reset ilen and imax for each row */ 649 for ( i=0; i<m; i++ ) { 650 ailen[i] = imax[i] = ai[i+1] - ai[i]; 651 } 652 a->nz = ai[m] + shift; 653 654 /* diagonals may have moved, so kill the diagonal pointers */ 655 if (fshift && a->diag) { 656 PetscFree(a->diag); 657 PLogObjectMemory(A,-(m+1)*sizeof(int)); 658 a->diag = 0; 659 } 660 PLogInfo(A,"MatAssemblyEnd_SeqAIJ:Matrix size: %d X %d; storage space: %d unneeded, %d used\n", 661 m,a->n,fshift,a->nz); 662 PLogInfo(A,"MatAssemblyEnd_SeqAIJ:Number of mallocs during MatSetValues is %d\n", 663 a->reallocs); 664 PLogInfo(A,"MatAssemblyEnd_SeqAIJ:Most nonzeros in any row is %d\n",rmax); 665 a->reallocs = 0; 666 A->info.nz_unneeded = (double)fshift; 667 668 /* check out for identical nodes. If found, use inode functions */ 669 ierr = Mat_AIJ_CheckInode(A); CHKERRQ(ierr); 670 return 0; 671 } 672 673 #undef __FUNC__ 674 #define __FUNC__ "MatZeroEntries_SeqAIJ" 675 int MatZeroEntries_SeqAIJ(Mat A) 676 { 677 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 678 PetscMemzero(a->a,(a->i[a->m]+a->indexshift)*sizeof(Scalar)); 679 return 0; 680 } 681 682 #undef __FUNC__ 683 #define __FUNC__ "MatDestroy_SeqAIJ" 684 int MatDestroy_SeqAIJ(PetscObject obj) 685 { 686 Mat A = (Mat) obj; 687 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 688 689 #if defined(PETSC_LOG) 690 PLogObjectState(obj,"Rows=%d, Cols=%d, NZ=%d",a->m,a->n,a->nz); 691 #endif 692 PetscFree(a->a); 693 if (!a->singlemalloc) { PetscFree(a->i); PetscFree(a->j);} 694 if (a->diag) PetscFree(a->diag); 695 if (a->ilen) PetscFree(a->ilen); 696 if (a->imax) PetscFree(a->imax); 697 if (a->solve_work) PetscFree(a->solve_work); 698 if (a->inode.size) PetscFree(a->inode.size); 699 PetscFree(a); 700 701 PLogObjectDestroy(A); 702 PetscHeaderDestroy(A); 703 return 0; 704 } 705 706 #undef __FUNC__ 707 #define __FUNC__ "MatCompress_SeqAIJ" 708 int MatCompress_SeqAIJ(Mat A) 709 { 710 return 0; 711 } 712 713 #undef __FUNC__ 714 #define __FUNC__ "MatSetOption_SeqAIJ" 715 int MatSetOption_SeqAIJ(Mat A,MatOption op) 716 { 717 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 718 if (op == MAT_ROW_ORIENTED) a->roworiented = 1; 719 else if (op == MAT_COLUMN_ORIENTED) a->roworiented = 0; 720 else if (op == MAT_COLUMNS_SORTED) a->sorted = 1; 721 else if (op == MAT_COLUMNS_UNSORTED) a->sorted = 0; 722 else if (op == MAT_NO_NEW_NONZERO_LOCATIONS) a->nonew = 1; 723 else if (op == MAT_NEW_NONZERO_LOCATION_ERROR) a->nonew = -1; 724 else if (op == MAT_NEW_NONZERO_ALLOCATION_ERROR) a->nonew = -2; 725 else if (op == MAT_YES_NEW_NONZERO_LOCATIONS) a->nonew = 0; 726 else if (op == MAT_ROWS_SORTED || 727 op == MAT_ROWS_UNSORTED || 728 op == MAT_SYMMETRIC || 729 op == MAT_STRUCTURALLY_SYMMETRIC || 730 op == MAT_YES_NEW_DIAGONALS || 731 op == MAT_IGNORE_OFF_PROC_ENTRIES) 732 PLogInfo(A,"Info:MatSetOption_SeqAIJ:Option ignored\n"); 733 else if (op == MAT_NO_NEW_DIAGONALS) 734 {SETERRQ(PETSC_ERR_SUP,0,"MAT_NO_NEW_DIAGONALS");} 735 else if (op == MAT_INODE_LIMIT_1) a->inode.limit = 1; 736 else if (op == MAT_INODE_LIMIT_2) a->inode.limit = 2; 737 else if (op == MAT_INODE_LIMIT_3) a->inode.limit = 3; 738 else if (op == MAT_INODE_LIMIT_4) a->inode.limit = 4; 739 else if (op == MAT_INODE_LIMIT_5) a->inode.limit = 5; 740 else 741 {SETERRQ(PETSC_ERR_SUP,0,"unknown option");} 742 return 0; 743 } 744 745 #undef __FUNC__ 746 #define __FUNC__ "MatGetDiagonal_SeqAIJ" 747 int MatGetDiagonal_SeqAIJ(Mat A,Vec v) 748 { 749 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 750 int i,j, n,shift = a->indexshift; 751 Scalar *x, zero = 0.0; 752 753 VecSet(&zero,v); 754 VecGetArray_Fast(v,x); VecGetLocalSize(v,&n); 755 if (n != a->m) SETERRQ(1,0,"Nonconforming matrix and vector"); 756 for ( i=0; i<a->m; i++ ) { 757 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 758 if (a->j[j]+shift == i) { 759 x[i] = a->a[j]; 760 break; 761 } 762 } 763 } 764 return 0; 765 } 766 767 /* -------------------------------------------------------*/ 768 /* Should check that shapes of vectors and matrices match */ 769 /* -------------------------------------------------------*/ 770 #undef __FUNC__ 771 #define __FUNC__ "MatMultTrans_SeqAIJ" 772 int MatMultTrans_SeqAIJ(Mat A,Vec xx,Vec yy) 773 { 774 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 775 Scalar *x, *y, *v, alpha; 776 int m = a->m, n, i, *idx, shift = a->indexshift; 777 778 VecGetArray_Fast(xx,x); VecGetArray_Fast(yy,y); 779 PetscMemzero(y,a->n*sizeof(Scalar)); 780 y = y + shift; /* shift for Fortran start by 1 indexing */ 781 for ( i=0; i<m; i++ ) { 782 idx = a->j + a->i[i] + shift; 783 v = a->a + a->i[i] + shift; 784 n = a->i[i+1] - a->i[i]; 785 alpha = x[i]; 786 while (n-->0) {y[*idx++] += alpha * *v++;} 787 } 788 PLogFlops(2*a->nz - a->n); 789 return 0; 790 } 791 792 #undef __FUNC__ 793 #define __FUNC__ "MatMultTransAdd_SeqAIJ" 794 int MatMultTransAdd_SeqAIJ(Mat A,Vec xx,Vec zz,Vec yy) 795 { 796 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 797 Scalar *x, *y, *v, alpha; 798 int m = a->m, n, i, *idx,shift = a->indexshift; 799 800 VecGetArray_Fast(xx,x); VecGetArray_Fast(yy,y); 801 if (zz != yy) VecCopy(zz,yy); 802 y = y + shift; /* shift for Fortran start by 1 indexing */ 803 for ( i=0; i<m; i++ ) { 804 idx = a->j + a->i[i] + shift; 805 v = a->a + a->i[i] + shift; 806 n = a->i[i+1] - a->i[i]; 807 alpha = x[i]; 808 while (n-->0) {y[*idx++] += alpha * *v++;} 809 } 810 PLogFlops(2*a->nz); 811 return 0; 812 } 813 814 #undef __FUNC__ 815 #define __FUNC__ "MatMult_SeqAIJ" 816 int MatMult_SeqAIJ(Mat A,Vec xx,Vec yy) 817 { 818 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 819 Scalar *x, *y, *v, sum; 820 int m = a->m, n, i, *idx, shift = a->indexshift,*ii,jrow,j; 821 822 VecGetArray_Fast(xx,x); VecGetArray_Fast(yy,y); 823 x = x + shift; /* shift for Fortran start by 1 indexing */ 824 idx = a->j; 825 v = a->a; 826 ii = a->i; 827 for ( i=0; i<m; i++ ) { 828 jrow = ii[i]; 829 n = ii[i+1] - jrow; 830 sum = 0.0; 831 /* while (n--) sum += *v++ * x[*idx++]; */ 832 for ( j=0; j<n; j++) { 833 sum += v[jrow]*x[idx[jrow]]; jrow++; 834 } 835 y[i] = sum; 836 } 837 PLogFlops(2*a->nz - m); 838 return 0; 839 } 840 841 #undef __FUNC__ 842 #define __FUNC__ "MatMultAdd_SeqAIJ" 843 int MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 844 { 845 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 846 Scalar *x, *y, *z, *v, sum; 847 int m = a->m, n, i, *idx, shift = a->indexshift,*ii,jrow,j; 848 849 850 VecGetArray_Fast(xx,x); VecGetArray_Fast(yy,y); VecGetArray_Fast(zz,z); 851 x = x + shift; /* shift for Fortran start by 1 indexing */ 852 idx = a->j; 853 v = a->a; 854 ii = a->i; 855 for ( i=0; i<m; i++ ) { 856 jrow = ii[i]; 857 n = ii[i+1] - jrow; 858 sum = y[i]; 859 /* while (n--) sum += *v++ * x[*idx++]; */ 860 for ( j=0; j<n; j++) { 861 sum += v[jrow]*x[idx[jrow]]; jrow++; 862 } 863 z[i] = sum; 864 } 865 PLogFlops(2*a->nz); 866 return 0; 867 } 868 869 /* 870 Adds diagonal pointers to sparse matrix structure. 871 */ 872 873 #undef __FUNC__ 874 #define __FUNC__ "MatMarkDiag_SeqAIJ" 875 int MatMarkDiag_SeqAIJ(Mat A) 876 { 877 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 878 int i,j, *diag, m = a->m,shift = a->indexshift; 879 880 diag = (int *) PetscMalloc( (m+1)*sizeof(int)); CHKPTRQ(diag); 881 PLogObjectMemory(A,(m+1)*sizeof(int)); 882 for ( i=0; i<a->m; i++ ) { 883 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 884 if (a->j[j]+shift == i) { 885 diag[i] = j - shift; 886 break; 887 } 888 } 889 } 890 a->diag = diag; 891 return 0; 892 } 893 894 #undef __FUNC__ 895 #define __FUNC__ "MatRelax_SeqAIJ" 896 int MatRelax_SeqAIJ(Mat A,Vec bb,double omega,MatSORType flag, 897 double fshift,int its,Vec xx) 898 { 899 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 900 Scalar *x, *b, *bs, d, *xs, sum, *v = a->a,*t,scale,*ts, *xb; 901 int ierr, *idx, *diag,n = a->n, m = a->m, i, shift = a->indexshift; 902 903 VecGetArray_Fast(xx,x); VecGetArray_Fast(bb,b); 904 if (!a->diag) {if ((ierr = MatMarkDiag_SeqAIJ(A))) return ierr;} 905 diag = a->diag; 906 xs = x + shift; /* shifted by one for index start of a or a->j*/ 907 if (flag == SOR_APPLY_UPPER) { 908 /* apply ( U + D/omega) to the vector */ 909 bs = b + shift; 910 for ( i=0; i<m; i++ ) { 911 d = fshift + a->a[diag[i] + shift]; 912 n = a->i[i+1] - diag[i] - 1; 913 idx = a->j + diag[i] + (!shift); 914 v = a->a + diag[i] + (!shift); 915 sum = b[i]*d/omega; 916 SPARSEDENSEDOT(sum,bs,v,idx,n); 917 x[i] = sum; 918 } 919 return 0; 920 } 921 if (flag == SOR_APPLY_LOWER) { 922 SETERRQ(1,0,"SOR_APPLY_LOWER is not done"); 923 } 924 else if (flag & SOR_EISENSTAT) { 925 /* Let A = L + U + D; where L is lower trianglar, 926 U is upper triangular, E is diagonal; This routine applies 927 928 (L + E)^{-1} A (U + E)^{-1} 929 930 to a vector efficiently using Eisenstat's trick. This is for 931 the case of SSOR preconditioner, so E is D/omega where omega 932 is the relaxation factor. 933 */ 934 t = (Scalar *) PetscMalloc( m*sizeof(Scalar) ); CHKPTRQ(t); 935 scale = (2.0/omega) - 1.0; 936 937 /* x = (E + U)^{-1} b */ 938 for ( i=m-1; i>=0; i-- ) { 939 d = fshift + a->a[diag[i] + shift]; 940 n = a->i[i+1] - diag[i] - 1; 941 idx = a->j + diag[i] + (!shift); 942 v = a->a + diag[i] + (!shift); 943 sum = b[i]; 944 SPARSEDENSEMDOT(sum,xs,v,idx,n); 945 x[i] = omega*(sum/d); 946 } 947 948 /* t = b - (2*E - D)x */ 949 v = a->a; 950 for ( i=0; i<m; i++ ) { t[i] = b[i] - scale*(v[*diag++ + shift])*x[i]; } 951 952 /* t = (E + L)^{-1}t */ 953 ts = t + shift; /* shifted by one for index start of a or a->j*/ 954 diag = a->diag; 955 for ( i=0; i<m; i++ ) { 956 d = fshift + a->a[diag[i]+shift]; 957 n = diag[i] - a->i[i]; 958 idx = a->j + a->i[i] + shift; 959 v = a->a + a->i[i] + shift; 960 sum = t[i]; 961 SPARSEDENSEMDOT(sum,ts,v,idx,n); 962 t[i] = omega*(sum/d); 963 } 964 965 /* x = x + t */ 966 for ( i=0; i<m; i++ ) { x[i] += t[i]; } 967 PetscFree(t); 968 return 0; 969 } 970 if (flag & SOR_ZERO_INITIAL_GUESS) { 971 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 972 for ( i=0; i<m; i++ ) { 973 d = fshift + a->a[diag[i]+shift]; 974 n = diag[i] - a->i[i]; 975 idx = a->j + a->i[i] + shift; 976 v = a->a + a->i[i] + shift; 977 sum = b[i]; 978 SPARSEDENSEMDOT(sum,xs,v,idx,n); 979 x[i] = omega*(sum/d); 980 } 981 xb = x; 982 } 983 else xb = b; 984 if ((flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) && 985 (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP)) { 986 for ( i=0; i<m; i++ ) { 987 x[i] *= a->a[diag[i]+shift]; 988 } 989 } 990 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 991 for ( i=m-1; i>=0; i-- ) { 992 d = fshift + a->a[diag[i] + shift]; 993 n = a->i[i+1] - diag[i] - 1; 994 idx = a->j + diag[i] + (!shift); 995 v = a->a + diag[i] + (!shift); 996 sum = xb[i]; 997 SPARSEDENSEMDOT(sum,xs,v,idx,n); 998 x[i] = omega*(sum/d); 999 } 1000 } 1001 its--; 1002 } 1003 while (its--) { 1004 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1005 for ( i=0; i<m; i++ ) { 1006 d = fshift + a->a[diag[i]+shift]; 1007 n = a->i[i+1] - a->i[i]; 1008 idx = a->j + a->i[i] + shift; 1009 v = a->a + a->i[i] + shift; 1010 sum = b[i]; 1011 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1012 x[i] = (1. - omega)*x[i] + omega*(sum + a->a[diag[i]+shift]*x[i])/d; 1013 } 1014 } 1015 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1016 for ( i=m-1; i>=0; i-- ) { 1017 d = fshift + a->a[diag[i] + shift]; 1018 n = a->i[i+1] - a->i[i]; 1019 idx = a->j + a->i[i] + shift; 1020 v = a->a + a->i[i] + shift; 1021 sum = b[i]; 1022 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1023 x[i] = (1. - omega)*x[i] + omega*(sum + a->a[diag[i]+shift]*x[i])/d; 1024 } 1025 } 1026 } 1027 return 0; 1028 } 1029 1030 #undef __FUNC__ 1031 #define __FUNC__ "MatGetInfo_SeqAIJ" 1032 int MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1033 { 1034 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1035 1036 info->rows_global = (double)a->m; 1037 info->columns_global = (double)a->n; 1038 info->rows_local = (double)a->m; 1039 info->columns_local = (double)a->n; 1040 info->block_size = 1.0; 1041 info->nz_allocated = (double)a->maxnz; 1042 info->nz_used = (double)a->nz; 1043 info->nz_unneeded = (double)(a->maxnz - a->nz); 1044 /* if (info->nz_unneeded != A->info.nz_unneeded) 1045 printf("space descrepancy: maxnz-nz = %d, nz_unneeded = %d\n",(int)info->nz_unneeded,(int)A->info.nz_unneeded); */ 1046 info->assemblies = (double)A->num_ass; 1047 info->mallocs = (double)a->reallocs; 1048 info->memory = A->mem; 1049 if (A->factor) { 1050 info->fill_ratio_given = A->info.fill_ratio_given; 1051 info->fill_ratio_needed = A->info.fill_ratio_needed; 1052 info->factor_mallocs = A->info.factor_mallocs; 1053 } else { 1054 info->fill_ratio_given = 0; 1055 info->fill_ratio_needed = 0; 1056 info->factor_mallocs = 0; 1057 } 1058 return 0; 1059 } 1060 1061 extern int MatLUFactorSymbolic_SeqAIJ(Mat,IS,IS,double,Mat*); 1062 extern int MatLUFactorNumeric_SeqAIJ(Mat,Mat*); 1063 extern int MatLUFactor_SeqAIJ(Mat,IS,IS,double); 1064 extern int MatSolve_SeqAIJ(Mat,Vec,Vec); 1065 extern int MatSolveAdd_SeqAIJ(Mat,Vec,Vec,Vec); 1066 extern int MatSolveTrans_SeqAIJ(Mat,Vec,Vec); 1067 extern int MatSolveTransAdd_SeqAIJ(Mat,Vec,Vec,Vec); 1068 1069 #undef __FUNC__ 1070 #define __FUNC__ "MatZeroRows_SeqAIJ" 1071 int MatZeroRows_SeqAIJ(Mat A,IS is,Scalar *diag) 1072 { 1073 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1074 int i,ierr,N, *rows,m = a->m - 1,shift = a->indexshift; 1075 1076 ierr = ISGetSize(is,&N); CHKERRQ(ierr); 1077 ierr = ISGetIndices(is,&rows); CHKERRQ(ierr); 1078 if (diag) { 1079 for ( i=0; i<N; i++ ) { 1080 if (rows[i] < 0 || rows[i] > m) SETERRQ(1,0,"row out of range"); 1081 if (a->ilen[rows[i]] > 0) { /* in case row was completely empty */ 1082 a->ilen[rows[i]] = 1; 1083 a->a[a->i[rows[i]]+shift] = *diag; 1084 a->j[a->i[rows[i]]+shift] = rows[i]+shift; 1085 } 1086 else { 1087 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],diag,INSERT_VALUES); 1088 CHKERRQ(ierr); 1089 } 1090 } 1091 } 1092 else { 1093 for ( i=0; i<N; i++ ) { 1094 if (rows[i] < 0 || rows[i] > m) SETERRQ(1,0,"row out of range"); 1095 a->ilen[rows[i]] = 0; 1096 } 1097 } 1098 ISRestoreIndices(is,&rows); 1099 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1100 return 0; 1101 } 1102 1103 #undef __FUNC__ 1104 #define __FUNC__ "MatGetSize_SeqAIJ" 1105 int MatGetSize_SeqAIJ(Mat A,int *m,int *n) 1106 { 1107 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1108 *m = a->m; *n = a->n; 1109 return 0; 1110 } 1111 1112 #undef __FUNC__ 1113 #define __FUNC__ "MatGetOwnershipRange_SeqAIJ" 1114 int MatGetOwnershipRange_SeqAIJ(Mat A,int *m,int *n) 1115 { 1116 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1117 *m = 0; *n = a->m; 1118 return 0; 1119 } 1120 1121 #undef __FUNC__ 1122 #define __FUNC__ "MatGetRow_SeqAIJ" 1123 int MatGetRow_SeqAIJ(Mat A,int row,int *nz,int **idx,Scalar **v) 1124 { 1125 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1126 int *itmp,i,shift = a->indexshift; 1127 1128 if (row < 0 || row >= a->m) SETERRQ(1,0,"Row out of range"); 1129 1130 *nz = a->i[row+1] - a->i[row]; 1131 if (v) *v = a->a + a->i[row] + shift; 1132 if (idx) { 1133 itmp = a->j + a->i[row] + shift; 1134 if (*nz && shift) { 1135 *idx = (int *) PetscMalloc( (*nz)*sizeof(int) ); CHKPTRQ(*idx); 1136 for ( i=0; i<(*nz); i++ ) {(*idx)[i] = itmp[i] + shift;} 1137 } else if (*nz) { 1138 *idx = itmp; 1139 } 1140 else *idx = 0; 1141 } 1142 return 0; 1143 } 1144 1145 #undef __FUNC__ 1146 #define __FUNC__ "MatRestoreRow_SeqAIJ" 1147 int MatRestoreRow_SeqAIJ(Mat A,int row,int *nz,int **idx,Scalar **v) 1148 { 1149 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1150 if (idx) {if (*idx && a->indexshift) PetscFree(*idx);} 1151 return 0; 1152 } 1153 1154 #undef __FUNC__ 1155 #define __FUNC__ "MatNorm_SeqAIJ" 1156 int MatNorm_SeqAIJ(Mat A,NormType type,double *norm) 1157 { 1158 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1159 Scalar *v = a->a; 1160 double sum = 0.0; 1161 int i, j,shift = a->indexshift; 1162 1163 if (type == NORM_FROBENIUS) { 1164 for (i=0; i<a->nz; i++ ) { 1165 #if defined(PETSC_COMPLEX) 1166 sum += real(conj(*v)*(*v)); v++; 1167 #else 1168 sum += (*v)*(*v); v++; 1169 #endif 1170 } 1171 *norm = sqrt(sum); 1172 } 1173 else if (type == NORM_1) { 1174 double *tmp; 1175 int *jj = a->j; 1176 tmp = (double *) PetscMalloc( (a->n+1)*sizeof(double) ); CHKPTRQ(tmp); 1177 PetscMemzero(tmp,a->n*sizeof(double)); 1178 *norm = 0.0; 1179 for ( j=0; j<a->nz; j++ ) { 1180 tmp[*jj++ + shift] += PetscAbsScalar(*v); v++; 1181 } 1182 for ( j=0; j<a->n; j++ ) { 1183 if (tmp[j] > *norm) *norm = tmp[j]; 1184 } 1185 PetscFree(tmp); 1186 } 1187 else if (type == NORM_INFINITY) { 1188 *norm = 0.0; 1189 for ( j=0; j<a->m; j++ ) { 1190 v = a->a + a->i[j] + shift; 1191 sum = 0.0; 1192 for ( i=0; i<a->i[j+1]-a->i[j]; i++ ) { 1193 sum += PetscAbsScalar(*v); v++; 1194 } 1195 if (sum > *norm) *norm = sum; 1196 } 1197 } 1198 else { 1199 SETERRQ(1,0,"No support for two norm yet"); 1200 } 1201 return 0; 1202 } 1203 1204 #undef __FUNC__ 1205 #define __FUNC__ "MatTranspose_SeqAIJ" 1206 int MatTranspose_SeqAIJ(Mat A,Mat *B) 1207 { 1208 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1209 Mat C; 1210 int i, ierr, *aj = a->j, *ai = a->i, m = a->m, len, *col; 1211 int shift = a->indexshift; 1212 Scalar *array = a->a; 1213 1214 if (B == PETSC_NULL && m != a->n) 1215 SETERRQ(1,0,"Square matrix only for in-place"); 1216 col = (int *) PetscMalloc((1+a->n)*sizeof(int)); CHKPTRQ(col); 1217 PetscMemzero(col,(1+a->n)*sizeof(int)); 1218 if (shift) { 1219 for ( i=0; i<ai[m]-1; i++ ) aj[i] -= 1; 1220 } 1221 for ( i=0; i<ai[m]+shift; i++ ) col[aj[i]] += 1; 1222 ierr = MatCreateSeqAIJ(A->comm,a->n,m,0,col,&C); CHKERRQ(ierr); 1223 PetscFree(col); 1224 for ( i=0; i<m; i++ ) { 1225 len = ai[i+1]-ai[i]; 1226 ierr = MatSetValues(C,len,aj,1,&i,array,INSERT_VALUES); CHKERRQ(ierr); 1227 array += len; aj += len; 1228 } 1229 if (shift) { 1230 for ( i=0; i<ai[m]-1; i++ ) aj[i] += 1; 1231 } 1232 1233 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1234 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1235 1236 if (B != PETSC_NULL) { 1237 *B = C; 1238 } else { 1239 /* This isn't really an in-place transpose */ 1240 PetscFree(a->a); 1241 if (!a->singlemalloc) {PetscFree(a->i); PetscFree(a->j);} 1242 if (a->diag) PetscFree(a->diag); 1243 if (a->ilen) PetscFree(a->ilen); 1244 if (a->imax) PetscFree(a->imax); 1245 if (a->solve_work) PetscFree(a->solve_work); 1246 if (a->inode.size) PetscFree(a->inode.size); 1247 PetscFree(a); 1248 PetscMemcpy(A,C,sizeof(struct _p_Mat)); 1249 PetscHeaderDestroy(C); 1250 } 1251 return 0; 1252 } 1253 1254 #undef __FUNC__ 1255 #define __FUNC__ "MatDiagonalScale_SeqAIJ" 1256 int MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 1257 { 1258 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1259 Scalar *l,*r,x,*v; 1260 int i,j,m = a->m, n = a->n, M, nz = a->nz, *jj,shift = a->indexshift; 1261 1262 if (ll) { 1263 /* The local size is used so that VecMPI can be passed to this routine 1264 by MatDiagonalScale_MPIAIJ */ 1265 VecGetLocalSize_Fast(ll,m); 1266 if (m != a->m) SETERRQ(1,0,"Left scaling vector wrong length"); 1267 VecGetArray_Fast(ll,l); 1268 v = a->a; 1269 for ( i=0; i<m; i++ ) { 1270 x = l[i]; 1271 M = a->i[i+1] - a->i[i]; 1272 for ( j=0; j<M; j++ ) { (*v++) *= x;} 1273 } 1274 PLogFlops(nz); 1275 } 1276 if (rr) { 1277 VecGetLocalSize_Fast(rr,n); 1278 if (n != a->n) SETERRQ(1,0,"Right scaling vector wrong length"); 1279 VecGetArray_Fast(rr,r); 1280 v = a->a; jj = a->j; 1281 for ( i=0; i<nz; i++ ) { 1282 (*v++) *= r[*jj++ + shift]; 1283 } 1284 PLogFlops(nz); 1285 } 1286 return 0; 1287 } 1288 1289 #undef __FUNC__ 1290 #define __FUNC__ "MatGetSubMatrix_SeqAIJ" 1291 int MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,MatGetSubMatrixCall scall,Mat *B) 1292 { 1293 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data,*c; 1294 int nznew, *smap, i, k, kstart, kend, ierr, oldcols = a->n,*lens; 1295 int row,mat_i,*mat_j,tcol,first,step,*mat_ilen; 1296 register int sum,lensi; 1297 int *irow, *icol, nrows, ncols, shift = a->indexshift,*ssmap; 1298 int *starts,*j_new,*i_new,*aj = a->j, *ai = a->i,ii,*ailen = a->ilen; 1299 Scalar *a_new,*mat_a; 1300 Mat C; 1301 1302 ierr = ISSorted(isrow,(PetscTruth*)&i); 1303 if (!i) SETERRQ(1,0,"ISrow is not sorted"); 1304 ierr = ISSorted(iscol,(PetscTruth*)&i); 1305 if (!i) SETERRQ(1,0,"IScol is not sorted"); 1306 1307 ierr = ISGetIndices(isrow,&irow); CHKERRQ(ierr); 1308 ierr = ISGetSize(isrow,&nrows); CHKERRQ(ierr); 1309 ierr = ISGetSize(iscol,&ncols); CHKERRQ(ierr); 1310 1311 if (ISStrideGetInfo(iscol,&first,&step) && step == 1) { /* no need to sort */ 1312 /* special case of contiguous rows */ 1313 lens = (int *) PetscMalloc((ncols+nrows+1)*sizeof(int)); CHKPTRQ(lens); 1314 starts = lens + ncols; 1315 /* loop over new rows determining lens and starting points */ 1316 for (i=0; i<nrows; i++) { 1317 kstart = ai[irow[i]]+shift; 1318 kend = kstart + ailen[irow[i]]; 1319 for ( k=kstart; k<kend; k++ ) { 1320 if (aj[k]+shift >= first) { 1321 starts[i] = k; 1322 break; 1323 } 1324 } 1325 sum = 0; 1326 while (k < kend) { 1327 if (aj[k++]+shift >= first+ncols) break; 1328 sum++; 1329 } 1330 lens[i] = sum; 1331 } 1332 /* create submatrix */ 1333 if (scall == MAT_REUSE_MATRIX) { 1334 int n_cols,n_rows; 1335 ierr = MatGetSize(*B,&n_rows,&n_cols); CHKERRQ(ierr); 1336 if (n_rows != nrows || n_cols != ncols) SETERRQ(1,0,""); 1337 ierr = MatZeroEntries(*B); CHKERRQ(ierr); 1338 C = *B; 1339 } 1340 else { 1341 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1342 } 1343 c = (Mat_SeqAIJ*) C->data; 1344 1345 /* loop over rows inserting into submatrix */ 1346 a_new = c->a; 1347 j_new = c->j; 1348 i_new = c->i; 1349 i_new[0] = -shift; 1350 for (i=0; i<nrows; i++) { 1351 ii = starts[i]; 1352 lensi = lens[i]; 1353 for ( k=0; k<lensi; k++ ) { 1354 *j_new++ = aj[ii+k] - first; 1355 } 1356 PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(Scalar)); 1357 a_new += lensi; 1358 i_new[i+1] = i_new[i] + lensi; 1359 c->ilen[i] = lensi; 1360 } 1361 PetscFree(lens); 1362 } 1363 else { 1364 ierr = ISGetIndices(iscol,&icol); CHKERRQ(ierr); 1365 smap = (int *) PetscMalloc((1+oldcols)*sizeof(int)); CHKPTRQ(smap); 1366 ssmap = smap + shift; 1367 lens = (int *) PetscMalloc((1+nrows)*sizeof(int)); CHKPTRQ(lens); 1368 PetscMemzero(smap,oldcols*sizeof(int)); 1369 for ( i=0; i<ncols; i++ ) smap[icol[i]] = i+1; 1370 /* determine lens of each row */ 1371 for (i=0; i<nrows; i++) { 1372 kstart = ai[irow[i]]+shift; 1373 kend = kstart + a->ilen[irow[i]]; 1374 lens[i] = 0; 1375 for ( k=kstart; k<kend; k++ ) { 1376 if (ssmap[aj[k]]) { 1377 lens[i]++; 1378 } 1379 } 1380 } 1381 /* Create and fill new matrix */ 1382 if (scall == MAT_REUSE_MATRIX) { 1383 c = (Mat_SeqAIJ *)((*B)->data); 1384 1385 if (c->m != nrows || c->n != ncols) SETERRQ(1,0,""); 1386 if (PetscMemcmp(c->ilen,lens, c->m *sizeof(int))) { 1387 SETERRQ(1,0,"Cannot reuse matrix. wrong no of nonzeros"); 1388 } 1389 PetscMemzero(c->ilen,c->m*sizeof(int)); 1390 C = *B; 1391 } 1392 else { 1393 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1394 } 1395 c = (Mat_SeqAIJ *)(C->data); 1396 for (i=0; i<nrows; i++) { 1397 row = irow[i]; 1398 nznew = 0; 1399 kstart = ai[row]+shift; 1400 kend = kstart + a->ilen[row]; 1401 mat_i = c->i[i]+shift; 1402 mat_j = c->j + mat_i; 1403 mat_a = c->a + mat_i; 1404 mat_ilen = c->ilen + i; 1405 for ( k=kstart; k<kend; k++ ) { 1406 if ((tcol=ssmap[a->j[k]])) { 1407 *mat_j++ = tcol - (!shift); 1408 *mat_a++ = a->a[k]; 1409 (*mat_ilen)++; 1410 1411 } 1412 } 1413 } 1414 /* Free work space */ 1415 ierr = ISRestoreIndices(iscol,&icol); CHKERRQ(ierr); 1416 PetscFree(smap); PetscFree(lens); 1417 } 1418 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1419 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1420 1421 ierr = ISRestoreIndices(isrow,&irow); CHKERRQ(ierr); 1422 *B = C; 1423 return 0; 1424 } 1425 1426 /* 1427 note: This can only work for identity for row and col. It would 1428 be good to check this and otherwise generate an error. 1429 */ 1430 #undef __FUNC__ 1431 #define __FUNC__ "MatILUFactor_SeqAIJ" 1432 int MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,double efill,int fill) 1433 { 1434 Mat_SeqAIJ *a = (Mat_SeqAIJ *) inA->data; 1435 int ierr; 1436 Mat outA; 1437 1438 if (fill != 0) SETERRQ(1,0,"Only fill=0 supported"); 1439 1440 outA = inA; 1441 inA->factor = FACTOR_LU; 1442 a->row = row; 1443 a->col = col; 1444 1445 if (!a->solve_work) { /* this matrix may have been factored before */ 1446 a->solve_work = (Scalar *) PetscMalloc( (a->m+1)*sizeof(Scalar)); CHKPTRQ(a->solve_work); 1447 } 1448 1449 if (!a->diag) { 1450 ierr = MatMarkDiag_SeqAIJ(inA); CHKERRQ(ierr); 1451 } 1452 ierr = MatLUFactorNumeric_SeqAIJ(inA,&outA); CHKERRQ(ierr); 1453 return 0; 1454 } 1455 1456 #include "pinclude/plapack.h" 1457 #undef __FUNC__ 1458 #define __FUNC__ "MatScale_SeqAIJ" 1459 int MatScale_SeqAIJ(Scalar *alpha,Mat inA) 1460 { 1461 Mat_SeqAIJ *a = (Mat_SeqAIJ *) inA->data; 1462 int one = 1; 1463 BLscal_( &a->nz, alpha, a->a, &one ); 1464 PLogFlops(a->nz); 1465 return 0; 1466 } 1467 1468 #undef __FUNC__ 1469 #define __FUNC__ "MatGetSubMatrices_SeqAIJ" 1470 int MatGetSubMatrices_SeqAIJ(Mat A,int n, IS *irow,IS *icol,MatGetSubMatrixCall scall, 1471 Mat **B) 1472 { 1473 int ierr,i; 1474 1475 if (scall == MAT_INITIAL_MATRIX) { 1476 *B = (Mat *) PetscMalloc( (n+1)*sizeof(Mat) ); CHKPTRQ(*B); 1477 } 1478 1479 for ( i=0; i<n; i++ ) { 1480 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],scall,&(*B)[i]);CHKERRQ(ierr); 1481 } 1482 return 0; 1483 } 1484 1485 #undef __FUNC__ 1486 #define __FUNC__ "MatGetBlockSize_SeqAIJ" 1487 int MatGetBlockSize_SeqAIJ(Mat A, int *bs) 1488 { 1489 *bs = 1; 1490 return 0; 1491 } 1492 1493 #undef __FUNC__ 1494 #define __FUNC__ "MatIncreaseOverlap_SeqAIJ" 1495 int MatIncreaseOverlap_SeqAIJ(Mat A, int is_max, IS *is, int ov) 1496 { 1497 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1498 int shift, row, i,j,k,l,m,n, *idx,ierr, *nidx, isz, val; 1499 int start, end, *ai, *aj; 1500 char *table; 1501 shift = a->indexshift; 1502 m = a->m; 1503 ai = a->i; 1504 aj = a->j+shift; 1505 1506 if (ov < 0) SETERRQ(1,0,"illegal overlap value used"); 1507 1508 table = (char *) PetscMalloc((m/BITSPERBYTE +1)*sizeof(char)); CHKPTRQ(table); 1509 nidx = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(nidx); 1510 1511 for ( i=0; i<is_max; i++ ) { 1512 /* Initialize the two local arrays */ 1513 isz = 0; 1514 PetscMemzero(table,(m/BITSPERBYTE +1)*sizeof(char)); 1515 1516 /* Extract the indices, assume there can be duplicate entries */ 1517 ierr = ISGetIndices(is[i],&idx); CHKERRQ(ierr); 1518 ierr = ISGetSize(is[i],&n); CHKERRQ(ierr); 1519 1520 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 1521 for ( j=0; j<n ; ++j){ 1522 if(!BT_LOOKUP(table, idx[j])) { nidx[isz++] = idx[j];} 1523 } 1524 ierr = ISRestoreIndices(is[i],&idx); CHKERRQ(ierr); 1525 ierr = ISDestroy(is[i]); CHKERRQ(ierr); 1526 1527 k = 0; 1528 for ( j=0; j<ov; j++){ /* for each overlap */ 1529 n = isz; 1530 for ( ; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 1531 row = nidx[k]; 1532 start = ai[row]; 1533 end = ai[row+1]; 1534 for ( l = start; l<end ; l++){ 1535 val = aj[l] + shift; 1536 if (!BT_LOOKUP(table,val)) {nidx[isz++] = val;} 1537 } 1538 } 1539 } 1540 ierr = ISCreateGeneral(PETSC_COMM_SELF, isz, nidx, (is+i)); CHKERRQ(ierr); 1541 } 1542 PetscFree(table); 1543 PetscFree(nidx); 1544 return 0; 1545 } 1546 1547 /* -------------------------------------------------------------- */ 1548 #undef __FUNC__ 1549 #define __FUNC__ "MatPermute_SeqAIJ" 1550 int MatPermute_SeqAIJ(Mat A, IS rowp, IS colp, Mat *B) 1551 { 1552 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1553 Scalar *vwork; 1554 int i, ierr, nz, m = a->m, n = a->n, *cwork; 1555 int *row,*col,*cnew,j,*lens; 1556 1557 ierr = ISGetIndices(rowp,&row); CHKERRQ(ierr); 1558 ierr = ISGetIndices(colp,&col); CHKERRQ(ierr); 1559 1560 /* determine lengths of permuted rows */ 1561 lens = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(lens); 1562 for (i=0; i<m; i++ ) { 1563 lens[row[i]] = a->i[i+1] - a->i[i]; 1564 } 1565 ierr = MatCreateSeqAIJ(A->comm,m,n,0,lens,B);CHKERRQ(ierr); 1566 PetscFree(lens); 1567 1568 cnew = (int *) PetscMalloc( n*sizeof(int) ); CHKPTRQ(cnew); 1569 for (i=0; i<m; i++) { 1570 ierr = MatGetRow(A,i,&nz,&cwork,&vwork); CHKERRQ(ierr); 1571 for (j=0; j<nz; j++ ) { cnew[j] = col[cwork[j]];} 1572 ierr = MatSetValues(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES); CHKERRQ(ierr); 1573 ierr = MatRestoreRow(A,i,&nz,&cwork,&vwork); CHKERRQ(ierr); 1574 } 1575 PetscFree(cnew); 1576 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1577 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1578 ierr = ISRestoreIndices(rowp,&row); CHKERRQ(ierr); 1579 ierr = ISRestoreIndices(colp,&col); CHKERRQ(ierr); 1580 return 0; 1581 } 1582 1583 #undef __FUNC__ 1584 #define __FUNC__ "MatPrintHelp_SeqAIJ" 1585 int MatPrintHelp_SeqAIJ(Mat A) 1586 { 1587 static int called = 0; 1588 MPI_Comm comm = A->comm; 1589 1590 if (called) return 0; else called = 1; 1591 PetscPrintf(comm," Options for MATSEQAIJ and MATMPIAIJ matrix formats (the defaults):\n"); 1592 PetscPrintf(comm," -mat_lu_pivotthreshold <threshold>: Set pivoting threshold\n"); 1593 PetscPrintf(comm," -mat_aij_oneindex: internal indices begin at 1 instead of the default 0.\n"); 1594 PetscPrintf(comm," -mat_aij_no_inode: Do not use inodes\n"); 1595 PetscPrintf(comm," -mat_aij_inode_limit <limit>: Set inode limit (max limit=5)\n"); 1596 #if defined(HAVE_ESSL) 1597 PetscPrintf(comm," -mat_aij_essl: Use IBM sparse LU factorization and solve.\n"); 1598 #endif 1599 return 0; 1600 } 1601 extern int MatEqual_SeqAIJ(Mat A,Mat B, PetscTruth* flg); 1602 extern int MatFDColoringCreate_SeqAIJ(Mat,ISColoring,MatFDColoring); 1603 extern int MatColoringPatch_SeqAIJ(Mat,int,int *,ISColoring *); 1604 1605 /* -------------------------------------------------------------------*/ 1606 static struct _MatOps MatOps = {MatSetValues_SeqAIJ, 1607 MatGetRow_SeqAIJ,MatRestoreRow_SeqAIJ, 1608 MatMult_SeqAIJ,MatMultAdd_SeqAIJ, 1609 MatMultTrans_SeqAIJ,MatMultTransAdd_SeqAIJ, 1610 MatSolve_SeqAIJ,MatSolveAdd_SeqAIJ, 1611 MatSolveTrans_SeqAIJ,MatSolveTransAdd_SeqAIJ, 1612 MatLUFactor_SeqAIJ,0, 1613 MatRelax_SeqAIJ, 1614 MatTranspose_SeqAIJ, 1615 MatGetInfo_SeqAIJ,MatEqual_SeqAIJ, 1616 MatGetDiagonal_SeqAIJ,MatDiagonalScale_SeqAIJ,MatNorm_SeqAIJ, 1617 0,MatAssemblyEnd_SeqAIJ, 1618 MatCompress_SeqAIJ, 1619 MatSetOption_SeqAIJ,MatZeroEntries_SeqAIJ,MatZeroRows_SeqAIJ, 1620 MatLUFactorSymbolic_SeqAIJ,MatLUFactorNumeric_SeqAIJ,0,0, 1621 MatGetSize_SeqAIJ,MatGetSize_SeqAIJ,MatGetOwnershipRange_SeqAIJ, 1622 MatILUFactorSymbolic_SeqAIJ,0, 1623 0,0, 1624 MatConvertSameType_SeqAIJ,0,0, 1625 MatILUFactor_SeqAIJ,0,0, 1626 MatGetSubMatrices_SeqAIJ,MatIncreaseOverlap_SeqAIJ, 1627 MatGetValues_SeqAIJ,0, 1628 MatPrintHelp_SeqAIJ, 1629 MatScale_SeqAIJ,0,0, 1630 MatILUDTFactor_SeqAIJ, 1631 MatGetBlockSize_SeqAIJ, 1632 MatGetRowIJ_SeqAIJ, 1633 MatRestoreRowIJ_SeqAIJ, 1634 MatGetColumnIJ_SeqAIJ, 1635 MatRestoreColumnIJ_SeqAIJ, 1636 MatFDColoringCreate_SeqAIJ, 1637 MatColoringPatch_SeqAIJ, 1638 0, 1639 MatPermute_SeqAIJ}; 1640 1641 extern int MatUseSuperLU_SeqAIJ(Mat); 1642 extern int MatUseEssl_SeqAIJ(Mat); 1643 extern int MatUseDXML_SeqAIJ(Mat); 1644 1645 #undef __FUNC__ 1646 #define __FUNC__ "MatCreateSeqAIJ" 1647 /*@C 1648 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 1649 (the default parallel PETSc format). For good matrix assembly performance 1650 the user should preallocate the matrix storage by setting the parameter nz 1651 (or the array nzz). By setting these parameters accurately, performance 1652 during matrix assembly can be increased by more than a factor of 50. 1653 1654 Input Parameters: 1655 . comm - MPI communicator, set to PETSC_COMM_SELF 1656 . m - number of rows 1657 . n - number of columns 1658 . nz - number of nonzeros per row (same for all rows) 1659 . nzz - array containing the number of nonzeros in the various rows 1660 (possibly different for each row) or PETSC_NULL 1661 1662 Output Parameter: 1663 . A - the matrix 1664 1665 Notes: 1666 The AIJ format (also called the Yale sparse matrix format or 1667 compressed row storage), is fully compatible with standard Fortran 77 1668 storage. That is, the stored row and column indices can begin at 1669 either one (as in Fortran) or zero. See the users' manual for details. 1670 1671 Specify the preallocated storage with either nz or nnz (not both). 1672 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1673 allocation. For large problems you MUST preallocate memory or you 1674 will get TERRIBLE performance, see the users' manual chapter on matrices. 1675 1676 By default, this format uses inodes (identical nodes) when possible, to 1677 improve numerical efficiency of Matrix vector products and solves. We 1678 search for consecutive rows with the same nonzero structure, thereby 1679 reusing matrix information to achieve increased efficiency. 1680 1681 Options Database Keys: 1682 $ -mat_aij_no_inode - Do not use inodes 1683 $ -mat_aij_inode_limit <limit> - Set inode limit. 1684 $ (max limit=5) 1685 $ -mat_aij_oneindex - Internally use indexing starting at 1 1686 $ rather than 0. Note: When calling MatSetValues(), 1687 $ the user still MUST index entries starting at 0! 1688 1689 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues() 1690 @*/ 1691 int MatCreateSeqAIJ(MPI_Comm comm,int m,int n,int nz,int *nnz, Mat *A) 1692 { 1693 Mat B; 1694 Mat_SeqAIJ *b; 1695 int i, len, ierr, flg,size; 1696 1697 MPI_Comm_size(comm,&size); 1698 if (size > 1) SETERRQ(1,0,"Comm must be of size 1"); 1699 1700 *A = 0; 1701 PetscHeaderCreate(B,_p_Mat,MAT_COOKIE,MATSEQAIJ,comm); 1702 PLogObjectCreate(B); 1703 B->data = (void *) (b = PetscNew(Mat_SeqAIJ)); CHKPTRQ(b); 1704 PetscMemzero(b,sizeof(Mat_SeqAIJ)); 1705 PetscMemcpy(&B->ops,&MatOps,sizeof(struct _MatOps)); 1706 B->destroy = MatDestroy_SeqAIJ; 1707 B->view = MatView_SeqAIJ; 1708 B->factor = 0; 1709 B->lupivotthreshold = 1.0; 1710 B->mapping = 0; 1711 ierr = OptionsGetDouble(PETSC_NULL,"-mat_lu_pivotthreshold",&B->lupivotthreshold, 1712 &flg); CHKERRQ(ierr); 1713 b->ilu_preserve_row_sums = PETSC_FALSE; 1714 ierr = OptionsHasName(PETSC_NULL,"-pc_ilu_preserve_row_sums", 1715 (int*) &b->ilu_preserve_row_sums); CHKERRQ(ierr); 1716 b->row = 0; 1717 b->col = 0; 1718 b->indexshift = 0; 1719 b->reallocs = 0; 1720 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_oneindex", &flg); CHKERRQ(ierr); 1721 if (flg) b->indexshift = -1; 1722 1723 b->m = m; B->m = m; B->M = m; 1724 b->n = n; B->n = n; B->N = n; 1725 b->imax = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(b->imax); 1726 if (nnz == PETSC_NULL) { 1727 if (nz == PETSC_DEFAULT) nz = 10; 1728 else if (nz <= 0) nz = 1; 1729 for ( i=0; i<m; i++ ) b->imax[i] = nz; 1730 nz = nz*m; 1731 } 1732 else { 1733 nz = 0; 1734 for ( i=0; i<m; i++ ) {b->imax[i] = nnz[i]; nz += nnz[i];} 1735 } 1736 1737 /* allocate the matrix space */ 1738 len = nz*(sizeof(int) + sizeof(Scalar)) + (b->m+1)*sizeof(int); 1739 b->a = (Scalar *) PetscMalloc( len ); CHKPTRQ(b->a); 1740 b->j = (int *) (b->a + nz); 1741 PetscMemzero(b->j,nz*sizeof(int)); 1742 b->i = b->j + nz; 1743 b->singlemalloc = 1; 1744 1745 b->i[0] = -b->indexshift; 1746 for (i=1; i<m+1; i++) { 1747 b->i[i] = b->i[i-1] + b->imax[i-1]; 1748 } 1749 1750 /* b->ilen will count nonzeros in each row so far. */ 1751 b->ilen = (int *) PetscMalloc((m+1)*sizeof(int)); 1752 PLogObjectMemory(B,len+2*(m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 1753 for ( i=0; i<b->m; i++ ) { b->ilen[i] = 0;} 1754 1755 b->nz = 0; 1756 b->maxnz = nz; 1757 b->sorted = 0; 1758 b->roworiented = 1; 1759 b->nonew = 0; 1760 b->diag = 0; 1761 b->solve_work = 0; 1762 b->spptr = 0; 1763 b->inode.node_count = 0; 1764 b->inode.size = 0; 1765 b->inode.limit = 5; 1766 b->inode.max_limit = 5; 1767 B->info.nz_unneeded = (double)b->maxnz; 1768 1769 *A = B; 1770 1771 /* SuperLU is not currently supported through PETSc */ 1772 #if defined(HAVE_SUPERLU) 1773 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_superlu", &flg); CHKERRQ(ierr); 1774 if (flg) { ierr = MatUseSuperLU_SeqAIJ(B); CHKERRQ(ierr); } 1775 #endif 1776 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_essl", &flg); CHKERRQ(ierr); 1777 if (flg) { ierr = MatUseEssl_SeqAIJ(B); CHKERRQ(ierr); } 1778 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_dxml", &flg); CHKERRQ(ierr); 1779 if (flg) { 1780 if (!b->indexshift) SETERRQ(1,0,"need -mat_aij_oneindex with -mat_aij_dxml"); 1781 ierr = MatUseDXML_SeqAIJ(B); CHKERRQ(ierr); 1782 } 1783 ierr = OptionsHasName(PETSC_NULL,"-help", &flg); CHKERRQ(ierr); 1784 if (flg) {ierr = MatPrintHelp(B); CHKERRQ(ierr); } 1785 return 0; 1786 } 1787 1788 #undef __FUNC__ 1789 #define __FUNC__ "MatConvertSameType_SeqAIJ" 1790 int MatConvertSameType_SeqAIJ(Mat A,Mat *B,int cpvalues) 1791 { 1792 Mat C; 1793 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ *) A->data; 1794 int i,len, m = a->m,shift = a->indexshift; 1795 1796 *B = 0; 1797 PetscHeaderCreate(C,_p_Mat,MAT_COOKIE,MATSEQAIJ,A->comm); 1798 PLogObjectCreate(C); 1799 C->data = (void *) (c = PetscNew(Mat_SeqAIJ)); CHKPTRQ(c); 1800 PetscMemcpy(&C->ops,&A->ops,sizeof(struct _MatOps)); 1801 C->destroy = MatDestroy_SeqAIJ; 1802 C->view = MatView_SeqAIJ; 1803 C->factor = A->factor; 1804 c->row = 0; 1805 c->col = 0; 1806 c->indexshift = shift; 1807 C->assembled = PETSC_TRUE; 1808 1809 c->m = C->m = a->m; 1810 c->n = C->n = a->n; 1811 C->M = a->m; 1812 C->N = a->n; 1813 1814 c->imax = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(c->imax); 1815 c->ilen = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(c->ilen); 1816 for ( i=0; i<m; i++ ) { 1817 c->imax[i] = a->imax[i]; 1818 c->ilen[i] = a->ilen[i]; 1819 } 1820 1821 /* allocate the matrix space */ 1822 c->singlemalloc = 1; 1823 len = (m+1)*sizeof(int)+(a->i[m])*(sizeof(Scalar)+sizeof(int)); 1824 c->a = (Scalar *) PetscMalloc( len ); CHKPTRQ(c->a); 1825 c->j = (int *) (c->a + a->i[m] + shift); 1826 c->i = c->j + a->i[m] + shift; 1827 PetscMemcpy(c->i,a->i,(m+1)*sizeof(int)); 1828 if (m > 0) { 1829 PetscMemcpy(c->j,a->j,(a->i[m]+shift)*sizeof(int)); 1830 if (cpvalues == COPY_VALUES) { 1831 PetscMemcpy(c->a,a->a,(a->i[m]+shift)*sizeof(Scalar)); 1832 } 1833 } 1834 1835 PLogObjectMemory(C,len+2*(m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 1836 c->sorted = a->sorted; 1837 c->roworiented = a->roworiented; 1838 c->nonew = a->nonew; 1839 c->ilu_preserve_row_sums = a->ilu_preserve_row_sums; 1840 1841 if (a->diag) { 1842 c->diag = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(c->diag); 1843 PLogObjectMemory(C,(m+1)*sizeof(int)); 1844 for ( i=0; i<m; i++ ) { 1845 c->diag[i] = a->diag[i]; 1846 } 1847 } 1848 else c->diag = 0; 1849 c->inode.limit = a->inode.limit; 1850 c->inode.max_limit = a->inode.max_limit; 1851 if (a->inode.size){ 1852 c->inode.size = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(c->inode.size); 1853 c->inode.node_count = a->inode.node_count; 1854 PetscMemcpy( c->inode.size, a->inode.size, (m+1)*sizeof(int)); 1855 } else { 1856 c->inode.size = 0; 1857 c->inode.node_count = 0; 1858 } 1859 c->nz = a->nz; 1860 c->maxnz = a->maxnz; 1861 c->solve_work = 0; 1862 c->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 1863 1864 *B = C; 1865 return 0; 1866 } 1867 1868 #undef __FUNC__ 1869 #define __FUNC__ "MatLoad_SeqAIJ" 1870 int MatLoad_SeqAIJ(Viewer viewer,MatType type,Mat *A) 1871 { 1872 Mat_SeqAIJ *a; 1873 Mat B; 1874 int i, nz, ierr, fd, header[4],size,*rowlengths = 0,M,N,shift; 1875 MPI_Comm comm; 1876 1877 PetscObjectGetComm((PetscObject) viewer,&comm); 1878 MPI_Comm_size(comm,&size); 1879 if (size > 1) SETERRQ(1,0,"view must have one processor"); 1880 ierr = ViewerBinaryGetDescriptor(viewer,&fd); CHKERRQ(ierr); 1881 ierr = PetscBinaryRead(fd,header,4,BINARY_INT); CHKERRQ(ierr); 1882 if (header[0] != MAT_COOKIE) SETERRQ(1,0,"not matrix object in file"); 1883 M = header[1]; N = header[2]; nz = header[3]; 1884 1885 /* read in row lengths */ 1886 rowlengths = (int*) PetscMalloc( M*sizeof(int) ); CHKPTRQ(rowlengths); 1887 ierr = PetscBinaryRead(fd,rowlengths,M,BINARY_INT); CHKERRQ(ierr); 1888 1889 /* create our matrix */ 1890 ierr = MatCreateSeqAIJ(comm,M,N,0,rowlengths,A); CHKERRQ(ierr); 1891 B = *A; 1892 a = (Mat_SeqAIJ *) B->data; 1893 shift = a->indexshift; 1894 1895 /* read in column indices and adjust for Fortran indexing*/ 1896 ierr = PetscBinaryRead(fd,a->j,nz,BINARY_INT); CHKERRQ(ierr); 1897 if (shift) { 1898 for ( i=0; i<nz; i++ ) { 1899 a->j[i] += 1; 1900 } 1901 } 1902 1903 /* read in nonzero values */ 1904 ierr = PetscBinaryRead(fd,a->a,nz,BINARY_SCALAR); CHKERRQ(ierr); 1905 1906 /* set matrix "i" values */ 1907 a->i[0] = -shift; 1908 for ( i=1; i<= M; i++ ) { 1909 a->i[i] = a->i[i-1] + rowlengths[i-1]; 1910 a->ilen[i-1] = rowlengths[i-1]; 1911 } 1912 PetscFree(rowlengths); 1913 1914 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1915 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1916 return 0; 1917 } 1918 1919 #undef __FUNC__ 1920 #define __FUNC__ "MatEqual_SeqAIJ" 1921 int MatEqual_SeqAIJ(Mat A,Mat B, PetscTruth* flg) 1922 { 1923 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)B->data; 1924 1925 if (B->type !=MATSEQAIJ)SETERRQ(1,0,"Matrices must be same type"); 1926 1927 /* If the matrix dimensions are not equal, or no of nonzeros or shift */ 1928 if ((a->m != b->m ) || (a->n !=b->n) ||( a->nz != b->nz)|| 1929 (a->indexshift != b->indexshift)) { 1930 *flg = PETSC_FALSE; return 0; 1931 } 1932 1933 /* if the a->i are the same */ 1934 if (PetscMemcmp(a->i,b->i, (a->n+1)*sizeof(int))) { 1935 *flg = PETSC_FALSE; return 0; 1936 } 1937 1938 /* if a->j are the same */ 1939 if (PetscMemcmp(a->j, b->j, (a->nz)*sizeof(int))) { 1940 *flg = PETSC_FALSE; return 0; 1941 } 1942 1943 /* if a->a are the same */ 1944 if (PetscMemcmp(a->a, b->a, (a->nz)*sizeof(Scalar))) { 1945 *flg = PETSC_FALSE; return 0; 1946 } 1947 *flg = PETSC_TRUE; 1948 return 0; 1949 1950 } 1951