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