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