1 #ifndef lint 2 static char vcid[] = "$Id: aij.c,v 1.220 1997/06/05 12:53:38 bsmith Exp curfman $"; 3 #endif 4 5 /* 6 Defines the basic matrix operations for the AIJ (compressed row) 7 matrix storage format. 8 */ 9 10 #include "pinclude/pviewer.h" 11 #include "sys.h" 12 #include "src/mat/impls/aij/seq/aij.h" 13 #include "src/vec/vecimpl.h" 14 #include "src/inline/spops.h" 15 #include "petsc.h" 16 #include "src/inline/bitarray.h" 17 18 /* 19 Basic AIJ format ILU based on drop tolerance 20 */ 21 #undef __FUNC__ 22 #define __FUNC__ "MatILUDTFactor_SeqAIJ" 23 int MatILUDTFactor_SeqAIJ(Mat A,double dt,int maxnz,IS row,IS col,Mat *fact) 24 { 25 /* Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; */ 26 int ierr = 1; 27 28 SETERRQ(ierr,0,"Not implemented"); 29 } 30 31 extern int MatToSymmetricIJ_SeqAIJ(int,int*,int*,int,int,int**,int**); 32 33 #undef __FUNC__ 34 #define __FUNC__ "MatGetRowIJ_SeqAIJ" 35 int MatGetRowIJ_SeqAIJ(Mat A,int oshift,PetscTruth symmetric,int *m,int **ia,int **ja, 36 PetscTruth *done) 37 { 38 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 39 int ierr,i,ishift; 40 41 *m = A->m; 42 if (!ia) return 0; 43 ishift = a->indexshift; 44 if (symmetric) { 45 ierr = MatToSymmetricIJ_SeqAIJ(a->m,a->i,a->j,ishift,oshift,ia,ja); CHKERRQ(ierr); 46 } else if (oshift == 0 && ishift == -1) { 47 int nz = a->i[a->m]; 48 /* malloc space and subtract 1 from i and j indices */ 49 *ia = (int *) PetscMalloc( (a->m+1)*sizeof(int) ); CHKPTRQ(*ia); 50 *ja = (int *) PetscMalloc( (nz+1)*sizeof(int) ); CHKPTRQ(*ja); 51 for ( i=0; i<nz; i++ ) (*ja)[i] = a->j[i] - 1; 52 for ( i=0; i<a->m+1; i++ ) (*ia)[i] = a->i[i] - 1; 53 } else if (oshift == 1 && ishift == 0) { 54 int nz = a->i[a->m] + 1; 55 /* malloc space and add 1 to i and j indices */ 56 *ia = (int *) PetscMalloc( (a->m+1)*sizeof(int) ); CHKPTRQ(*ia); 57 *ja = (int *) PetscMalloc( (nz+1)*sizeof(int) ); CHKPTRQ(*ja); 58 for ( i=0; i<nz; i++ ) (*ja)[i] = a->j[i] + 1; 59 for ( i=0; i<a->m+1; i++ ) (*ia)[i] = a->i[i] + 1; 60 } else { 61 *ia = a->i; *ja = a->j; 62 } 63 64 return 0; 65 } 66 67 #undef __FUNC__ 68 #define __FUNC__ "MatRestoreRowIJ_SeqAIJ" 69 int MatRestoreRowIJ_SeqAIJ(Mat A,int oshift,PetscTruth symmetric,int *n,int **ia,int **ja, 70 PetscTruth *done) 71 { 72 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 73 int ishift = a->indexshift; 74 75 if (!ia) return 0; 76 if (symmetric || (oshift == 0 && ishift == -1) || (oshift == 1 && ishift == 0)) { 77 PetscFree(*ia); 78 PetscFree(*ja); 79 } 80 return 0; 81 } 82 83 #undef __FUNC__ 84 #define __FUNC__ "MatGetColumnIJ_SeqAIJ" 85 int MatGetColumnIJ_SeqAIJ(Mat A,int oshift,PetscTruth symmetric,int *nn,int **ia,int **ja, 86 PetscTruth *done) 87 { 88 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 89 int ierr,i,ishift = a->indexshift,*collengths,*cia,*cja,n = A->n,m = A->m; 90 int nz = a->i[m]+ishift,row,*jj,mr,col; 91 92 *nn = A->n; 93 if (!ia) return 0; 94 if (symmetric) { 95 ierr = MatToSymmetricIJ_SeqAIJ(a->m,a->i,a->j,ishift,oshift,ia,ja); CHKERRQ(ierr); 96 } else { 97 collengths = (int *) PetscMalloc( (n+1)*sizeof(int) ); CHKPTRQ(collengths); 98 PetscMemzero(collengths,n*sizeof(int)); 99 cia = (int *) PetscMalloc( (n+1)*sizeof(int) ); CHKPTRQ(cia); 100 cja = (int *) PetscMalloc( (nz+1)*sizeof(int) ); CHKPTRQ(cja); 101 jj = a->j; 102 for ( i=0; i<nz; i++ ) { 103 collengths[jj[i] + ishift]++; 104 } 105 cia[0] = oshift; 106 for ( i=0; i<n; i++) { 107 cia[i+1] = cia[i] + collengths[i]; 108 } 109 PetscMemzero(collengths,n*sizeof(int)); 110 jj = a->j; 111 for ( row=0; row<m; row++ ) { 112 mr = a->i[row+1] - a->i[row]; 113 for ( i=0; i<mr; i++ ) { 114 col = *jj++ + ishift; 115 cja[cia[col] + collengths[col]++ - oshift] = row + oshift; 116 } 117 } 118 PetscFree(collengths); 119 *ia = cia; *ja = cja; 120 } 121 122 return 0; 123 } 124 125 #undef __FUNC__ 126 #define __FUNC__ "MatRestoreColumnIJ_SeqAIJ" 127 int MatRestoreColumnIJ_SeqAIJ(Mat A,int oshift,PetscTruth symmetric,int *n,int **ia, 128 int **ja,PetscTruth *done) 129 { 130 if (!ia) return 0; 131 132 PetscFree(*ia); 133 PetscFree(*ja); 134 135 return 0; 136 } 137 138 #define CHUNKSIZE 15 139 140 #undef __FUNC__ 141 #define __FUNC__ "MatSetValues_SeqAIJ" 142 int MatSetValues_SeqAIJ(Mat A,int m,int *im,int n,int *in,Scalar *v,InsertMode is) 143 { 144 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 145 int *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax, N, sorted = a->sorted; 146 int *imax = a->imax, *ai = a->i, *ailen = a->ilen,roworiented = a->roworiented; 147 int *aj = a->j, nonew = a->nonew,shift = a->indexshift; 148 Scalar *ap,value, *aa = a->a; 149 150 for ( k=0; k<m; k++ ) { /* loop over added rows */ 151 row = im[k]; 152 #if defined(PETSC_BOPT_g) 153 if (row < 0) SETERRQ(1,0,"Negative row"); 154 if (row >= a->m) SETERRQ(1,0,"Row too large"); 155 #endif 156 rp = aj + ai[row] + shift; ap = aa + ai[row] + shift; 157 rmax = imax[row]; nrow = ailen[row]; 158 low = 0; 159 for ( l=0; l<n; l++ ) { /* loop over added columns */ 160 #if defined(PETSC_BOPT_g) 161 if (in[l] < 0) SETERRQ(1,0,"Negative column"); 162 if (in[l] >= a->n) SETERRQ(1,0,"Column too large"); 163 #endif 164 col = in[l] - shift; 165 if (roworiented) { 166 value = *v++; 167 } 168 else { 169 value = v[k + l*m]; 170 } 171 if (!sorted) low = 0; high = nrow; 172 while (high-low > 5) { 173 t = (low+high)/2; 174 if (rp[t] > col) high = t; 175 else low = t; 176 } 177 for ( i=low; i<high; i++ ) { 178 if (rp[i] > col) break; 179 if (rp[i] == col) { 180 if (is == ADD_VALUES) ap[i] += value; 181 else ap[i] = value; 182 goto noinsert; 183 } 184 } 185 if (nonew == 1) goto noinsert; 186 else if (nonew == -1) SETERRQ(1,0,"Inserting a new nonzero in the matrix"); 187 if (nrow >= rmax) { 188 /* there is no extra room in row, therefore enlarge */ 189 int new_nz = ai[a->m] + CHUNKSIZE,len,*new_i,*new_j; 190 Scalar *new_a; 191 192 if (nonew == -2) SETERRQ(1,0,"Inserting a new nonzero in the matrix"); 193 194 /* malloc new storage space */ 195 len = new_nz*(sizeof(int)+sizeof(Scalar))+(a->m+1)*sizeof(int); 196 new_a = (Scalar *) PetscMalloc( len ); CHKPTRQ(new_a); 197 new_j = (int *) (new_a + new_nz); 198 new_i = new_j + new_nz; 199 200 /* copy over old data into new slots */ 201 for ( ii=0; ii<row+1; ii++ ) {new_i[ii] = ai[ii];} 202 for ( ii=row+1; ii<a->m+1; ii++ ) {new_i[ii] = ai[ii]+CHUNKSIZE;} 203 PetscMemcpy(new_j,aj,(ai[row]+nrow+shift)*sizeof(int)); 204 len = (new_nz - CHUNKSIZE - ai[row] - nrow - shift); 205 PetscMemcpy(new_j+ai[row]+shift+nrow+CHUNKSIZE,aj+ai[row]+shift+nrow, 206 len*sizeof(int)); 207 PetscMemcpy(new_a,aa,(ai[row]+nrow+shift)*sizeof(Scalar)); 208 PetscMemcpy(new_a+ai[row]+shift+nrow+CHUNKSIZE,aa+ai[row]+shift+nrow, 209 len*sizeof(Scalar)); 210 /* free up old matrix storage */ 211 PetscFree(a->a); 212 if (!a->singlemalloc) {PetscFree(a->i);PetscFree(a->j);} 213 aa = a->a = new_a; ai = a->i = new_i; aj = a->j = new_j; 214 a->singlemalloc = 1; 215 216 rp = aj + ai[row] + shift; ap = aa + ai[row] + shift; 217 rmax = imax[row] = imax[row] + CHUNKSIZE; 218 PLogObjectMemory(A,CHUNKSIZE*(sizeof(int) + sizeof(Scalar))); 219 a->maxnz += CHUNKSIZE; 220 a->reallocs++; 221 } 222 N = nrow++ - 1; a->nz++; 223 /* shift up all the later entries in this row */ 224 for ( ii=N; ii>=i; ii-- ) { 225 rp[ii+1] = rp[ii]; 226 ap[ii+1] = ap[ii]; 227 } 228 rp[i] = col; 229 ap[i] = value; 230 noinsert:; 231 low = i + 1; 232 } 233 ailen[row] = nrow; 234 } 235 return 0; 236 } 237 238 #undef __FUNC__ 239 #define __FUNC__ "MatGetValues_SeqAIJ" 240 int MatGetValues_SeqAIJ(Mat A,int m,int *im,int n,int *in,Scalar *v) 241 { 242 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 243 int *rp, k, low, high, t, row, nrow, i, col, l, *aj = a->j; 244 int *ai = a->i, *ailen = a->ilen, shift = a->indexshift; 245 Scalar *ap, *aa = a->a, zero = 0.0; 246 247 for ( k=0; k<m; k++ ) { /* loop over rows */ 248 row = im[k]; 249 if (row < 0) SETERRQ(1,0,"Negative row"); 250 if (row >= a->m) SETERRQ(1,0,"Row too large"); 251 rp = aj + ai[row] + shift; ap = aa + ai[row] + shift; 252 nrow = ailen[row]; 253 for ( l=0; l<n; l++ ) { /* loop over columns */ 254 if (in[l] < 0) SETERRQ(1,0,"Negative column"); 255 if (in[l] >= a->n) SETERRQ(1,0,"Column too large"); 256 col = in[l] - shift; 257 high = nrow; low = 0; /* assume unsorted */ 258 while (high-low > 5) { 259 t = (low+high)/2; 260 if (rp[t] > col) high = t; 261 else low = t; 262 } 263 for ( i=low; i<high; i++ ) { 264 if (rp[i] > col) break; 265 if (rp[i] == col) { 266 *v++ = ap[i]; 267 goto finished; 268 } 269 } 270 *v++ = zero; 271 finished:; 272 } 273 } 274 return 0; 275 } 276 277 278 #undef __FUNC__ 279 #define __FUNC__ "MatView_SeqAIJ_Binary" 280 extern int MatView_SeqAIJ_Binary(Mat A,Viewer viewer) 281 { 282 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 283 int i, fd, *col_lens, ierr; 284 285 ierr = ViewerBinaryGetDescriptor(viewer,&fd); CHKERRQ(ierr); 286 col_lens = (int *) PetscMalloc( (4+a->m)*sizeof(int) ); CHKPTRQ(col_lens); 287 col_lens[0] = MAT_COOKIE; 288 col_lens[1] = a->m; 289 col_lens[2] = a->n; 290 col_lens[3] = a->nz; 291 292 /* store lengths of each row and write (including header) to file */ 293 for ( i=0; i<a->m; i++ ) { 294 col_lens[4+i] = a->i[i+1] - a->i[i]; 295 } 296 ierr = PetscBinaryWrite(fd,col_lens,4+a->m,BINARY_INT,1); CHKERRQ(ierr); 297 PetscFree(col_lens); 298 299 /* store column indices (zero start index) */ 300 if (a->indexshift) { 301 for ( i=0; i<a->nz; i++ ) a->j[i]--; 302 } 303 ierr = PetscBinaryWrite(fd,a->j,a->nz,BINARY_INT,0); CHKERRQ(ierr); 304 if (a->indexshift) { 305 for ( i=0; i<a->nz; i++ ) a->j[i]++; 306 } 307 308 /* store nonzero values */ 309 ierr = PetscBinaryWrite(fd,a->a,a->nz,BINARY_SCALAR,0); CHKERRQ(ierr); 310 return 0; 311 } 312 313 #undef __FUNC__ 314 #define __FUNC__ "MatView_SeqAIJ_ASCII" 315 extern int MatView_SeqAIJ_ASCII(Mat A,Viewer viewer) 316 { 317 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 318 int ierr, i,j, m = a->m, shift = a->indexshift, format, flg1,flg2; 319 FILE *fd; 320 char *outputname; 321 322 ierr = ViewerASCIIGetPointer(viewer,&fd); CHKERRQ(ierr); 323 ierr = ViewerFileGetOutputname_Private(viewer,&outputname); CHKERRQ(ierr); 324 ierr = ViewerGetFormat(viewer,&format); 325 if (format == VIEWER_FORMAT_ASCII_INFO) { 326 return 0; 327 } 328 else if (format == VIEWER_FORMAT_ASCII_INFO_LONG) { 329 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_no_inode",&flg1); CHKERRQ(ierr); 330 ierr = OptionsHasName(PETSC_NULL,"-mat_no_unroll",&flg2); CHKERRQ(ierr); 331 if (flg1 || flg2) fprintf(fd," not using I-node routines\n"); 332 else fprintf(fd," using I-node routines: found %d nodes, limit used is %d\n", 333 a->inode.node_count,a->inode.limit); 334 } 335 else if (format == VIEWER_FORMAT_ASCII_MATLAB) { 336 fprintf(fd,"%% Size = %d %d \n",m,a->n); 337 fprintf(fd,"%% Nonzeros = %d \n",a->nz); 338 fprintf(fd,"zzz = zeros(%d,3);\n",a->nz); 339 fprintf(fd,"zzz = [\n"); 340 341 for (i=0; i<m; i++) { 342 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 343 #if defined(PETSC_COMPLEX) 344 fprintf(fd,"%d %d %18.16e + %18.16e i \n",i+1,a->j[j]+!shift,real(a->a[j]), 345 imag(a->a[j])); 346 #else 347 fprintf(fd,"%d %d %18.16e\n", i+1, a->j[j]+!shift, a->a[j]); 348 #endif 349 } 350 } 351 fprintf(fd,"];\n %s = spconvert(zzz);\n",outputname); 352 } 353 else if (format == VIEWER_FORMAT_ASCII_COMMON) { 354 for ( i=0; i<m; i++ ) { 355 fprintf(fd,"row %d:",i); 356 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 357 #if defined(PETSC_COMPLEX) 358 if (imag(a->a[j]) != 0.0 && real(a->a[j]) != 0.0) 359 fprintf(fd," %d %g + %g i",a->j[j]+shift,real(a->a[j]),imag(a->a[j])); 360 else if (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 extern 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 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 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 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 643 #if defined(PETSC_LOG) 644 PLogObjectState(obj,"Rows=%d, Cols=%d, NZ=%d",a->m,a->n,a->nz); 645 #endif 646 PetscFree(a->a); 647 if (!a->singlemalloc) { PetscFree(a->i); PetscFree(a->j);} 648 if (a->diag) PetscFree(a->diag); 649 if (a->ilen) PetscFree(a->ilen); 650 if (a->imax) PetscFree(a->imax); 651 if (a->solve_work) PetscFree(a->solve_work); 652 if (a->inode.size) PetscFree(a->inode.size); 653 PetscFree(a); 654 655 PLogObjectDestroy(A); 656 PetscHeaderDestroy(A); 657 return 0; 658 } 659 660 #undef __FUNC__ 661 #define __FUNC__ "MatCompress_SeqAIJ" 662 int MatCompress_SeqAIJ(Mat A) 663 { 664 return 0; 665 } 666 667 #undef __FUNC__ 668 #define __FUNC__ "MatSetOption_SeqAIJ" 669 int MatSetOption_SeqAIJ(Mat A,MatOption op) 670 { 671 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 672 if (op == MAT_ROW_ORIENTED) a->roworiented = 1; 673 else if (op == MAT_COLUMN_ORIENTED) a->roworiented = 0; 674 else if (op == MAT_COLUMNS_SORTED) a->sorted = 1; 675 else if (op == MAT_COLUMNS_UNSORTED) a->sorted = 0; 676 else if (op == MAT_NO_NEW_NONZERO_LOCATIONS) a->nonew = 1; 677 else if (op == MAT_NEW_NONZERO_LOCATION_ERROR) a->nonew = -1; 678 else if (op == MAT_NEW_NONZERO_ALLOCATION_ERROR) a->nonew = -2; 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 _p_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,_p_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 _p_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,_p_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 _p_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