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