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