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