1 2 3 #ifndef lint 4 static char vcid[] = "$Id: aij.c,v 1.197 1996/11/29 22:20:00 curfman Exp curfman $"; 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 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 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_COLUMNS_UNSORTED) a->sorted = 0; 598 else if (op == MAT_NO_NEW_NONZERO_LOCATIONS) a->nonew = 1; 599 else if (op == MAT_YES_NEW_NONZERO_LOCATIONS) a->nonew = 0; 600 else if (op == MAT_ROWS_SORTED || 601 op == MAT_ROWS_UNSORTED || 602 op == MAT_SYMMETRIC || 603 op == MAT_STRUCTURALLY_SYMMETRIC || 604 op == MAT_YES_NEW_DIAGONALS || 605 op == MAT_IGNORE_OFF_PROCESSOR_ENTRIES) 606 PLogInfo(A,"Info:MatSetOption_SeqAIJ:Option ignored\n"); 607 else if (op == MAT_NO_NEW_DIAGONALS) 608 {SETERRQ(PETSC_ERR_SUP,"MatSetOption_SeqAIJ:MAT_NO_NEW_DIAGONALS");} 609 else if (op == MAT_INODE_LIMIT_1) a->inode.limit = 1; 610 else if (op == MAT_INODE_LIMIT_2) a->inode.limit = 2; 611 else if (op == MAT_INODE_LIMIT_3) a->inode.limit = 3; 612 else if (op == MAT_INODE_LIMIT_4) a->inode.limit = 4; 613 else if (op == MAT_INODE_LIMIT_5) a->inode.limit = 5; 614 else 615 {SETERRQ(PETSC_ERR_SUP,"MatSetOption_SeqAIJ:unknown option");} 616 return 0; 617 } 618 619 static int MatGetDiagonal_SeqAIJ(Mat A,Vec v) 620 { 621 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 622 int i,j, n,shift = a->indexshift; 623 Scalar *x, zero = 0.0; 624 625 VecSet(&zero,v); 626 VecGetArray_Fast(v,x); VecGetLocalSize(v,&n); 627 if (n != a->m) SETERRQ(1,"MatGetDiagonal_SeqAIJ:Nonconforming matrix and vector"); 628 for ( i=0; i<a->m; i++ ) { 629 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 630 if (a->j[j]+shift == i) { 631 x[i] = a->a[j]; 632 break; 633 } 634 } 635 } 636 return 0; 637 } 638 639 /* -------------------------------------------------------*/ 640 /* Should check that shapes of vectors and matrices match */ 641 /* -------------------------------------------------------*/ 642 int MatMultTrans_SeqAIJ(Mat A,Vec xx,Vec yy) 643 { 644 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 645 Scalar *x, *y, *v, alpha; 646 int m = a->m, n, i, *idx, shift = a->indexshift; 647 648 VecGetArray_Fast(xx,x); VecGetArray_Fast(yy,y); 649 PetscMemzero(y,a->n*sizeof(Scalar)); 650 y = y + shift; /* shift for Fortran start by 1 indexing */ 651 for ( i=0; i<m; i++ ) { 652 idx = a->j + a->i[i] + shift; 653 v = a->a + a->i[i] + shift; 654 n = a->i[i+1] - a->i[i]; 655 alpha = x[i]; 656 while (n-->0) {y[*idx++] += alpha * *v++;} 657 } 658 PLogFlops(2*a->nz - a->n); 659 return 0; 660 } 661 662 int MatMultTransAdd_SeqAIJ(Mat A,Vec xx,Vec zz,Vec yy) 663 { 664 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 665 Scalar *x, *y, *v, alpha; 666 int m = a->m, n, i, *idx,shift = a->indexshift; 667 668 VecGetArray_Fast(xx,x); VecGetArray_Fast(yy,y); 669 if (zz != yy) VecCopy(zz,yy); 670 y = y + shift; /* shift for Fortran start by 1 indexing */ 671 for ( i=0; i<m; i++ ) { 672 idx = a->j + a->i[i] + shift; 673 v = a->a + a->i[i] + shift; 674 n = a->i[i+1] - a->i[i]; 675 alpha = x[i]; 676 while (n-->0) {y[*idx++] += alpha * *v++;} 677 } 678 PLogFlops(2*a->nz); 679 return 0; 680 } 681 682 int MatMult_SeqAIJ(Mat A,Vec xx,Vec yy) 683 { 684 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 685 Scalar *x, *y, *v, sum; 686 int m = a->m, n, i, *idx, shift = a->indexshift,*ii; 687 688 VecGetArray_Fast(xx,x); VecGetArray_Fast(yy,y); 689 x = x + shift; /* shift for Fortran start by 1 indexing */ 690 idx = a->j; 691 v = a->a; 692 ii = a->i; 693 for ( i=0; i<m; i++ ) { 694 n = ii[1] - ii[0]; ii++; 695 sum = 0.0; 696 /* SPARSEDENSEDOT(sum,x,v,idx,n); */ 697 /* for ( j=n-1; j>-1; j--) sum += v[j]*x[idx[j]]; */ 698 while (n--) sum += *v++ * x[*idx++]; 699 y[i] = sum; 700 } 701 PLogFlops(2*a->nz - m); 702 return 0; 703 } 704 705 int MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 706 { 707 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 708 Scalar *x, *y, *z, *v, sum; 709 int m = a->m, n, i, *idx, shift = a->indexshift,*ii; 710 711 VecGetArray_Fast(xx,x); VecGetArray_Fast(yy,y); VecGetArray_Fast(zz,z); 712 x = x + shift; /* shift for Fortran start by 1 indexing */ 713 idx = a->j; 714 v = a->a; 715 ii = a->i; 716 for ( i=0; i<m; i++ ) { 717 n = ii[1] - ii[0]; ii++; 718 sum = y[i]; 719 /* SPARSEDENSEDOT(sum,x,v,idx,n); */ 720 while (n--) sum += *v++ * x[*idx++]; 721 z[i] = sum; 722 } 723 PLogFlops(2*a->nz); 724 return 0; 725 } 726 727 /* 728 Adds diagonal pointers to sparse matrix structure. 729 */ 730 731 int MatMarkDiag_SeqAIJ(Mat A) 732 { 733 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 734 int i,j, *diag, m = a->m,shift = a->indexshift; 735 736 diag = (int *) PetscMalloc( (m+1)*sizeof(int)); CHKPTRQ(diag); 737 PLogObjectMemory(A,(m+1)*sizeof(int)); 738 for ( i=0; i<a->m; i++ ) { 739 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 740 if (a->j[j]+shift == i) { 741 diag[i] = j - shift; 742 break; 743 } 744 } 745 } 746 a->diag = diag; 747 return 0; 748 } 749 750 int MatRelax_SeqAIJ(Mat A,Vec bb,double omega,MatSORType flag, 751 double fshift,int its,Vec xx) 752 { 753 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 754 Scalar *x, *b, *bs, d, *xs, sum, *v = a->a,*t,scale,*ts, *xb; 755 int ierr, *idx, *diag,n = a->n, m = a->m, i, shift = a->indexshift; 756 757 VecGetArray_Fast(xx,x); VecGetArray_Fast(bb,b); 758 if (!a->diag) {if ((ierr = MatMarkDiag_SeqAIJ(A))) return ierr;} 759 diag = a->diag; 760 xs = x + shift; /* shifted by one for index start of a or a->j*/ 761 if (flag == SOR_APPLY_UPPER) { 762 /* apply ( U + D/omega) to the vector */ 763 bs = b + shift; 764 for ( i=0; i<m; i++ ) { 765 d = fshift + a->a[diag[i] + shift]; 766 n = a->i[i+1] - diag[i] - 1; 767 idx = a->j + diag[i] + (!shift); 768 v = a->a + diag[i] + (!shift); 769 sum = b[i]*d/omega; 770 SPARSEDENSEDOT(sum,bs,v,idx,n); 771 x[i] = sum; 772 } 773 return 0; 774 } 775 if (flag == SOR_APPLY_LOWER) { 776 SETERRQ(1,"MatRelax_SeqAIJ:SOR_APPLY_LOWER is not done"); 777 } 778 else if (flag & SOR_EISENSTAT) { 779 /* Let A = L + U + D; where L is lower trianglar, 780 U is upper triangular, E is diagonal; This routine applies 781 782 (L + E)^{-1} A (U + E)^{-1} 783 784 to a vector efficiently using Eisenstat's trick. This is for 785 the case of SSOR preconditioner, so E is D/omega where omega 786 is the relaxation factor. 787 */ 788 t = (Scalar *) PetscMalloc( m*sizeof(Scalar) ); CHKPTRQ(t); 789 scale = (2.0/omega) - 1.0; 790 791 /* x = (E + U)^{-1} b */ 792 for ( i=m-1; i>=0; i-- ) { 793 d = fshift + a->a[diag[i] + shift]; 794 n = a->i[i+1] - diag[i] - 1; 795 idx = a->j + diag[i] + (!shift); 796 v = a->a + diag[i] + (!shift); 797 sum = b[i]; 798 SPARSEDENSEMDOT(sum,xs,v,idx,n); 799 x[i] = omega*(sum/d); 800 } 801 802 /* t = b - (2*E - D)x */ 803 v = a->a; 804 for ( i=0; i<m; i++ ) { t[i] = b[i] - scale*(v[*diag++ + shift])*x[i]; } 805 806 /* t = (E + L)^{-1}t */ 807 ts = t + shift; /* shifted by one for index start of a or a->j*/ 808 diag = a->diag; 809 for ( i=0; i<m; i++ ) { 810 d = fshift + a->a[diag[i]+shift]; 811 n = diag[i] - a->i[i]; 812 idx = a->j + a->i[i] + shift; 813 v = a->a + a->i[i] + shift; 814 sum = t[i]; 815 SPARSEDENSEMDOT(sum,ts,v,idx,n); 816 t[i] = omega*(sum/d); 817 } 818 819 /* x = x + t */ 820 for ( i=0; i<m; i++ ) { x[i] += t[i]; } 821 PetscFree(t); 822 return 0; 823 } 824 if (flag & SOR_ZERO_INITIAL_GUESS) { 825 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 826 for ( i=0; i<m; i++ ) { 827 d = fshift + a->a[diag[i]+shift]; 828 n = diag[i] - a->i[i]; 829 idx = a->j + a->i[i] + shift; 830 v = a->a + a->i[i] + shift; 831 sum = b[i]; 832 SPARSEDENSEMDOT(sum,xs,v,idx,n); 833 x[i] = omega*(sum/d); 834 } 835 xb = x; 836 } 837 else xb = b; 838 if ((flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) && 839 (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP)) { 840 for ( i=0; i<m; i++ ) { 841 x[i] *= a->a[diag[i]+shift]; 842 } 843 } 844 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 845 for ( i=m-1; i>=0; i-- ) { 846 d = fshift + a->a[diag[i] + shift]; 847 n = a->i[i+1] - diag[i] - 1; 848 idx = a->j + diag[i] + (!shift); 849 v = a->a + diag[i] + (!shift); 850 sum = xb[i]; 851 SPARSEDENSEMDOT(sum,xs,v,idx,n); 852 x[i] = omega*(sum/d); 853 } 854 } 855 its--; 856 } 857 while (its--) { 858 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 859 for ( i=0; i<m; i++ ) { 860 d = fshift + a->a[diag[i]+shift]; 861 n = a->i[i+1] - a->i[i]; 862 idx = a->j + a->i[i] + shift; 863 v = a->a + a->i[i] + shift; 864 sum = b[i]; 865 SPARSEDENSEMDOT(sum,xs,v,idx,n); 866 x[i] = (1. - omega)*x[i] + omega*(sum + a->a[diag[i]+shift]*x[i])/d; 867 } 868 } 869 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 870 for ( i=m-1; i>=0; i-- ) { 871 d = fshift + a->a[diag[i] + shift]; 872 n = a->i[i+1] - a->i[i]; 873 idx = a->j + a->i[i] + shift; 874 v = a->a + a->i[i] + shift; 875 sum = b[i]; 876 SPARSEDENSEMDOT(sum,xs,v,idx,n); 877 x[i] = (1. - omega)*x[i] + omega*(sum + a->a[diag[i]+shift]*x[i])/d; 878 } 879 } 880 } 881 return 0; 882 } 883 884 static int MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 885 { 886 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 887 888 info->rows_global = (double)a->m; 889 info->columns_global = (double)a->n; 890 info->rows_local = (double)a->m; 891 info->columns_local = (double)a->n; 892 info->block_size = 1.0; 893 info->nz_allocated = (double)a->maxnz; 894 info->nz_used = (double)a->nz; 895 info->nz_unneeded = (double)(a->maxnz - a->nz); 896 /* if (info->nz_unneeded != A->info.nz_unneeded) 897 printf("space descrepancy: maxnz-nz = %d, nz_unneeded = %d\n",(int)info->nz_unneeded,(int)A->info.nz_unneeded); */ 898 info->assemblies = (double)A->num_ass; 899 info->mallocs = (double)a->reallocs; 900 info->memory = A->mem; 901 if (A->factor) { 902 info->fill_ratio_given = A->info.fill_ratio_given; 903 info->fill_ratio_needed = A->info.fill_ratio_needed; 904 info->factor_mallocs = A->info.factor_mallocs; 905 } else { 906 info->fill_ratio_given = 0; 907 info->fill_ratio_needed = 0; 908 info->factor_mallocs = 0; 909 } 910 return 0; 911 } 912 913 extern int MatLUFactorSymbolic_SeqAIJ(Mat,IS,IS,double,Mat*); 914 extern int MatLUFactorNumeric_SeqAIJ(Mat,Mat*); 915 extern int MatLUFactor_SeqAIJ(Mat,IS,IS,double); 916 extern int MatSolve_SeqAIJ(Mat,Vec,Vec); 917 extern int MatSolveAdd_SeqAIJ(Mat,Vec,Vec,Vec); 918 extern int MatSolveTrans_SeqAIJ(Mat,Vec,Vec); 919 extern int MatSolveTransAdd_SeqAIJ(Mat,Vec,Vec,Vec); 920 921 static int MatZeroRows_SeqAIJ(Mat A,IS is,Scalar *diag) 922 { 923 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 924 int i,ierr,N, *rows,m = a->m - 1,shift = a->indexshift; 925 926 ierr = ISGetSize(is,&N); CHKERRQ(ierr); 927 ierr = ISGetIndices(is,&rows); CHKERRQ(ierr); 928 if (diag) { 929 for ( i=0; i<N; i++ ) { 930 if (rows[i] < 0 || rows[i] > m) SETERRQ(1,"MatZeroRows_SeqAIJ:row out of range"); 931 if (a->ilen[rows[i]] > 0) { /* in case row was completely empty */ 932 a->ilen[rows[i]] = 1; 933 a->a[a->i[rows[i]]+shift] = *diag; 934 a->j[a->i[rows[i]]+shift] = rows[i]+shift; 935 } 936 else { 937 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],diag,INSERT_VALUES); 938 CHKERRQ(ierr); 939 } 940 } 941 } 942 else { 943 for ( i=0; i<N; i++ ) { 944 if (rows[i] < 0 || rows[i] > m) SETERRQ(1,"MatZeroRows_SeqAIJ:row out of range"); 945 a->ilen[rows[i]] = 0; 946 } 947 } 948 ISRestoreIndices(is,&rows); 949 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 950 return 0; 951 } 952 953 static int MatGetSize_SeqAIJ(Mat A,int *m,int *n) 954 { 955 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 956 *m = a->m; *n = a->n; 957 return 0; 958 } 959 960 static int MatGetOwnershipRange_SeqAIJ(Mat A,int *m,int *n) 961 { 962 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 963 *m = 0; *n = a->m; 964 return 0; 965 } 966 int MatGetRow_SeqAIJ(Mat A,int row,int *nz,int **idx,Scalar **v) 967 { 968 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 969 int *itmp,i,shift = a->indexshift; 970 971 if (row < 0 || row >= a->m) SETERRQ(1,"MatGetRow_SeqAIJ:Row out of range"); 972 973 *nz = a->i[row+1] - a->i[row]; 974 if (v) *v = a->a + a->i[row] + shift; 975 if (idx) { 976 itmp = a->j + a->i[row] + shift; 977 if (*nz && shift) { 978 *idx = (int *) PetscMalloc( (*nz)*sizeof(int) ); CHKPTRQ(*idx); 979 for ( i=0; i<(*nz); i++ ) {(*idx)[i] = itmp[i] + shift;} 980 } else if (*nz) { 981 *idx = itmp; 982 } 983 else *idx = 0; 984 } 985 return 0; 986 } 987 988 int MatRestoreRow_SeqAIJ(Mat A,int row,int *nz,int **idx,Scalar **v) 989 { 990 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 991 if (idx) {if (*idx && a->indexshift) PetscFree(*idx);} 992 return 0; 993 } 994 995 static int MatNorm_SeqAIJ(Mat A,NormType type,double *norm) 996 { 997 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 998 Scalar *v = a->a; 999 double sum = 0.0; 1000 int i, j,shift = a->indexshift; 1001 1002 if (type == NORM_FROBENIUS) { 1003 for (i=0; i<a->nz; i++ ) { 1004 #if defined(PETSC_COMPLEX) 1005 sum += real(conj(*v)*(*v)); v++; 1006 #else 1007 sum += (*v)*(*v); v++; 1008 #endif 1009 } 1010 *norm = sqrt(sum); 1011 } 1012 else if (type == NORM_1) { 1013 double *tmp; 1014 int *jj = a->j; 1015 tmp = (double *) PetscMalloc( a->n*sizeof(double) ); CHKPTRQ(tmp); 1016 PetscMemzero(tmp,a->n*sizeof(double)); 1017 *norm = 0.0; 1018 for ( j=0; j<a->nz; j++ ) { 1019 tmp[*jj++ + shift] += PetscAbsScalar(*v); v++; 1020 } 1021 for ( j=0; j<a->n; j++ ) { 1022 if (tmp[j] > *norm) *norm = tmp[j]; 1023 } 1024 PetscFree(tmp); 1025 } 1026 else if (type == NORM_INFINITY) { 1027 *norm = 0.0; 1028 for ( j=0; j<a->m; j++ ) { 1029 v = a->a + a->i[j] + shift; 1030 sum = 0.0; 1031 for ( i=0; i<a->i[j+1]-a->i[j]; i++ ) { 1032 sum += PetscAbsScalar(*v); v++; 1033 } 1034 if (sum > *norm) *norm = sum; 1035 } 1036 } 1037 else { 1038 SETERRQ(1,"MatNorm_SeqAIJ:No support for two norm yet"); 1039 } 1040 return 0; 1041 } 1042 1043 static int MatTranspose_SeqAIJ(Mat A,Mat *B) 1044 { 1045 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1046 Mat C; 1047 int i, ierr, *aj = a->j, *ai = a->i, m = a->m, len, *col; 1048 int shift = a->indexshift; 1049 Scalar *array = a->a; 1050 1051 if (B == PETSC_NULL && m != a->n) 1052 SETERRQ(1,"MatTranspose_SeqAIJ:Square matrix only for in-place"); 1053 col = (int *) PetscMalloc((1+a->n)*sizeof(int)); CHKPTRQ(col); 1054 PetscMemzero(col,(1+a->n)*sizeof(int)); 1055 if (shift) { 1056 for ( i=0; i<ai[m]-1; i++ ) aj[i] -= 1; 1057 } 1058 for ( i=0; i<ai[m]+shift; i++ ) col[aj[i]] += 1; 1059 ierr = MatCreateSeqAIJ(A->comm,a->n,m,0,col,&C); CHKERRQ(ierr); 1060 PetscFree(col); 1061 for ( i=0; i<m; i++ ) { 1062 len = ai[i+1]-ai[i]; 1063 ierr = MatSetValues(C,len,aj,1,&i,array,INSERT_VALUES); CHKERRQ(ierr); 1064 array += len; aj += len; 1065 } 1066 if (shift) { 1067 for ( i=0; i<ai[m]-1; i++ ) aj[i] += 1; 1068 } 1069 1070 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1071 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1072 1073 if (B != PETSC_NULL) { 1074 *B = C; 1075 } else { 1076 /* This isn't really an in-place transpose */ 1077 PetscFree(a->a); 1078 if (!a->singlemalloc) {PetscFree(a->i); PetscFree(a->j);} 1079 if (a->diag) PetscFree(a->diag); 1080 if (a->ilen) PetscFree(a->ilen); 1081 if (a->imax) PetscFree(a->imax); 1082 if (a->solve_work) PetscFree(a->solve_work); 1083 if (a->inode.size) PetscFree(a->inode.size); 1084 PetscFree(a); 1085 PetscMemcpy(A,C,sizeof(struct _Mat)); 1086 PetscHeaderDestroy(C); 1087 } 1088 return 0; 1089 } 1090 1091 static int MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 1092 { 1093 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1094 Scalar *l,*r,x,*v; 1095 int i,j,m = a->m, n = a->n, M, nz = a->nz, *jj,shift = a->indexshift; 1096 1097 if (ll) { 1098 /* The local size is used so that VecMPI can be passed to this routine 1099 by MatDiagonalScale_MPIAIJ */ 1100 VecGetLocalSize_Fast(ll,m); 1101 if (m != a->m) SETERRQ(1,"MatDiagonalScale_SeqAIJ:Left scaling vector wrong length"); 1102 VecGetArray_Fast(ll,l); 1103 v = a->a; 1104 for ( i=0; i<m; i++ ) { 1105 x = l[i]; 1106 M = a->i[i+1] - a->i[i]; 1107 for ( j=0; j<M; j++ ) { (*v++) *= x;} 1108 } 1109 PLogFlops(nz); 1110 } 1111 if (rr) { 1112 VecGetLocalSize_Fast(rr,n); 1113 if (n != a->n) SETERRQ(1,"MatDiagonalScale_SeqAIJ:Right scaling vector wrong length"); 1114 VecGetArray_Fast(rr,r); 1115 v = a->a; jj = a->j; 1116 for ( i=0; i<nz; i++ ) { 1117 (*v++) *= r[*jj++ + shift]; 1118 } 1119 PLogFlops(nz); 1120 } 1121 return 0; 1122 } 1123 1124 static int MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,MatGetSubMatrixCall scall,Mat *B) 1125 { 1126 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data,*c; 1127 int nznew, *smap, i, k, kstart, kend, ierr, oldcols = a->n,*lens; 1128 int row,mat_i,*mat_j,tcol,first,step,*mat_ilen; 1129 register int sum,lensi; 1130 int *irow, *icol, nrows, ncols, shift = a->indexshift,*ssmap; 1131 int *starts,*j_new,*i_new,*aj = a->j, *ai = a->i,ii,*ailen = a->ilen; 1132 Scalar *a_new,*mat_a; 1133 Mat C; 1134 1135 ierr = ISSorted(isrow,(PetscTruth*)&i); 1136 if (!i) SETERRQ(1,"MatGetSubmatrices_SeqAIJ:ISrow is not sorted"); 1137 ierr = ISSorted(iscol,(PetscTruth*)&i); 1138 if (!i) SETERRQ(1,"MatGetSubmatrices_SeqAIJ:IScol is not sorted"); 1139 1140 ierr = ISGetIndices(isrow,&irow); CHKERRQ(ierr); 1141 ierr = ISGetSize(isrow,&nrows); CHKERRQ(ierr); 1142 ierr = ISGetSize(iscol,&ncols); CHKERRQ(ierr); 1143 1144 if (ISStrideGetInfo(iscol,&first,&step) && step == 1) { /* no need to sort */ 1145 /* special case of contiguous rows */ 1146 lens = (int *) PetscMalloc((ncols+nrows+1)*sizeof(int)); CHKPTRQ(lens); 1147 starts = lens + ncols; 1148 /* loop over new rows determining lens and starting points */ 1149 for (i=0; i<nrows; i++) { 1150 kstart = ai[irow[i]]+shift; 1151 kend = kstart + ailen[irow[i]]; 1152 for ( k=kstart; k<kend; k++ ) { 1153 if (aj[k]+shift >= first) { 1154 starts[i] = k; 1155 break; 1156 } 1157 } 1158 sum = 0; 1159 while (k < kend) { 1160 if (aj[k++]+shift >= first+ncols) break; 1161 sum++; 1162 } 1163 lens[i] = sum; 1164 } 1165 /* create submatrix */ 1166 if (scall == MAT_REUSE_MATRIX) { 1167 int n_cols,n_rows; 1168 ierr = MatGetSize(*B,&n_rows,&n_cols); CHKERRQ(ierr); 1169 if (n_rows != nrows || n_cols != ncols) SETERRQ(1,"MatGetSubMatrix_SeqAIJ:"); 1170 ierr = MatZeroEntries(*B); CHKERRQ(ierr); 1171 C = *B; 1172 } 1173 else { 1174 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1175 } 1176 c = (Mat_SeqAIJ*) C->data; 1177 1178 /* loop over rows inserting into submatrix */ 1179 a_new = c->a; 1180 j_new = c->j; 1181 i_new = c->i; 1182 i_new[0] = -shift; 1183 for (i=0; i<nrows; i++) { 1184 ii = starts[i]; 1185 lensi = lens[i]; 1186 for ( k=0; k<lensi; k++ ) { 1187 *j_new++ = aj[ii+k] - first; 1188 } 1189 PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(Scalar)); 1190 a_new += lensi; 1191 i_new[i+1] = i_new[i] + lensi; 1192 c->ilen[i] = lensi; 1193 } 1194 PetscFree(lens); 1195 } 1196 else { 1197 ierr = ISGetIndices(iscol,&icol); CHKERRQ(ierr); 1198 smap = (int *) PetscMalloc((1+oldcols)*sizeof(int)); CHKPTRQ(smap); 1199 ssmap = smap + shift; 1200 lens = (int *) PetscMalloc((1+nrows)*sizeof(int)); CHKPTRQ(lens); 1201 PetscMemzero(smap,oldcols*sizeof(int)); 1202 for ( i=0; i<ncols; i++ ) smap[icol[i]] = i+1; 1203 /* determine lens of each row */ 1204 for (i=0; i<nrows; i++) { 1205 kstart = ai[irow[i]]+shift; 1206 kend = kstart + a->ilen[irow[i]]; 1207 lens[i] = 0; 1208 for ( k=kstart; k<kend; k++ ) { 1209 if (ssmap[aj[k]]) { 1210 lens[i]++; 1211 } 1212 } 1213 } 1214 /* Create and fill new matrix */ 1215 if (scall == MAT_REUSE_MATRIX) { 1216 c = (Mat_SeqAIJ *)((*B)->data); 1217 1218 if (c->m != nrows || c->n != ncols) SETERRQ(1,"MatGetSubMatrix_SeqAIJ:"); 1219 if (PetscMemcmp(c->ilen,lens, c->m *sizeof(int))) { 1220 SETERRQ(1,"MatGetSubmatrices_SeqAIJ:Cannot reuse matrix. wrong no of nonzeros"); 1221 } 1222 PetscMemzero(c->ilen,c->m*sizeof(int)); 1223 C = *B; 1224 } 1225 else { 1226 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1227 } 1228 c = (Mat_SeqAIJ *)(C->data); 1229 for (i=0; i<nrows; i++) { 1230 row = irow[i]; 1231 nznew = 0; 1232 kstart = ai[row]+shift; 1233 kend = kstart + a->ilen[row]; 1234 mat_i = c->i[i]+shift; 1235 mat_j = c->j + mat_i; 1236 mat_a = c->a + mat_i; 1237 mat_ilen = c->ilen + i; 1238 for ( k=kstart; k<kend; k++ ) { 1239 if ((tcol=ssmap[a->j[k]])) { 1240 *mat_j++ = tcol - (!shift); 1241 *mat_a++ = a->a[k]; 1242 (*mat_ilen)++; 1243 1244 } 1245 } 1246 } 1247 /* Free work space */ 1248 ierr = ISRestoreIndices(iscol,&icol); CHKERRQ(ierr); 1249 PetscFree(smap); PetscFree(lens); 1250 } 1251 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1252 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1253 1254 ierr = ISRestoreIndices(isrow,&irow); CHKERRQ(ierr); 1255 *B = C; 1256 return 0; 1257 } 1258 1259 /* 1260 note: This can only work for identity for row and col. It would 1261 be good to check this and otherwise generate an error. 1262 */ 1263 static int MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,double efill,int fill) 1264 { 1265 Mat_SeqAIJ *a = (Mat_SeqAIJ *) inA->data; 1266 int ierr; 1267 Mat outA; 1268 1269 if (fill != 0) SETERRQ(1,"MatILUFactor_SeqAIJ:Only fill=0 supported"); 1270 1271 outA = inA; 1272 inA->factor = FACTOR_LU; 1273 a->row = row; 1274 a->col = col; 1275 1276 a->solve_work = (Scalar *) PetscMalloc( (a->m+1)*sizeof(Scalar)); CHKPTRQ(a->solve_work); 1277 1278 if (!a->diag) { 1279 ierr = MatMarkDiag_SeqAIJ(inA); CHKERRQ(ierr); 1280 } 1281 ierr = MatLUFactorNumeric_SeqAIJ(inA,&outA); CHKERRQ(ierr); 1282 return 0; 1283 } 1284 1285 #include "pinclude/plapack.h" 1286 static int MatScale_SeqAIJ(Scalar *alpha,Mat inA) 1287 { 1288 Mat_SeqAIJ *a = (Mat_SeqAIJ *) inA->data; 1289 int one = 1; 1290 BLscal_( &a->nz, alpha, a->a, &one ); 1291 PLogFlops(a->nz); 1292 return 0; 1293 } 1294 1295 static int MatGetSubMatrices_SeqAIJ(Mat A,int n, IS *irow,IS *icol,MatGetSubMatrixCall scall, 1296 Mat **B) 1297 { 1298 int ierr,i; 1299 1300 if (scall == MAT_INITIAL_MATRIX) { 1301 *B = (Mat *) PetscMalloc( (n+1)*sizeof(Mat) ); CHKPTRQ(*B); 1302 } 1303 1304 for ( i=0; i<n; i++ ) { 1305 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],scall,&(*B)[i]);CHKERRQ(ierr); 1306 } 1307 return 0; 1308 } 1309 1310 static int MatGetBlockSize_SeqAIJ(Mat A, int *bs) 1311 { 1312 *bs = 1; 1313 return 0; 1314 } 1315 1316 static int MatIncreaseOverlap_SeqAIJ(Mat A, int is_max, IS *is, int ov) 1317 { 1318 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1319 int shift, row, i,j,k,l,m,n, *idx,ierr, *nidx, isz, val; 1320 int start, end, *ai, *aj; 1321 char *table; 1322 shift = a->indexshift; 1323 m = a->m; 1324 ai = a->i; 1325 aj = a->j+shift; 1326 1327 if (ov < 0) SETERRQ(1,"MatIncreaseOverlap_SeqAIJ: illegal overlap value used"); 1328 1329 table = (char *) PetscMalloc((m/BITSPERBYTE +1)*sizeof(char)); CHKPTRQ(table); 1330 nidx = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(nidx); 1331 1332 for ( i=0; i<is_max; i++ ) { 1333 /* Initialise the two local arrays */ 1334 isz = 0; 1335 PetscMemzero(table,(m/BITSPERBYTE +1)*sizeof(char)); 1336 1337 /* Extract the indices, assume there can be duplicate entries */ 1338 ierr = ISGetIndices(is[i],&idx); CHKERRQ(ierr); 1339 ierr = ISGetSize(is[i],&n); CHKERRQ(ierr); 1340 1341 /* Enter these into the temp arrays i.e mark table[row], enter row into new index */ 1342 for ( j=0; j<n ; ++j){ 1343 if(!BT_LOOKUP(table, idx[j])) { nidx[isz++] = idx[j];} 1344 } 1345 ierr = ISRestoreIndices(is[i],&idx); CHKERRQ(ierr); 1346 ierr = ISDestroy(is[i]); CHKERRQ(ierr); 1347 1348 k = 0; 1349 for ( j=0; j<ov; j++){ /* for each overlap*/ 1350 n = isz; 1351 for ( ; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 1352 row = nidx[k]; 1353 start = ai[row]; 1354 end = ai[row+1]; 1355 for ( l = start; l<end ; l++){ 1356 val = aj[l] + shift; 1357 if (!BT_LOOKUP(table,val)) {nidx[isz++] = val;} 1358 } 1359 } 1360 } 1361 ierr = ISCreateGeneral(MPI_COMM_SELF, isz, nidx, (is+i)); CHKERRQ(ierr); 1362 } 1363 PetscFree(table); 1364 PetscFree(nidx); 1365 return 0; 1366 } 1367 1368 int MatPrintHelp_SeqAIJ(Mat A) 1369 { 1370 static int called = 0; 1371 MPI_Comm comm = A->comm; 1372 1373 if (called) return 0; else called = 1; 1374 PetscPrintf(comm," Options for MATSEQAIJ and MATMPIAIJ matrix formats (the defaults):\n"); 1375 PetscPrintf(comm," -mat_lu_pivotthreshold <threshold>\n"); 1376 PetscPrintf(comm," -mat_aij_oneindex - internal indices begin at 1 instead of the default 0.\n"); 1377 PetscPrintf(comm," -mat_aij_no_inode - Do not use inodes\n"); 1378 PetscPrintf(comm," -mat_aij_inode_limit <limit> - Set inode limit (max limit=5)\n"); 1379 #if defined(HAVE_ESSL) 1380 PetscPrintf(comm," -mat_aij_essl - Use IBM sparse LU factorization and solve.\n"); 1381 #endif 1382 return 0; 1383 } 1384 static int MatEqual_SeqAIJ(Mat A,Mat B, PetscTruth* flg); 1385 extern int MatFDColoringCreate_SeqAIJ(Mat,ISColoring,MatFDColoring); 1386 extern int MatColoringPatch_SeqAIJ(Mat,int,int *,ISColoring *); 1387 1388 /* -------------------------------------------------------------------*/ 1389 static struct _MatOps MatOps = {MatSetValues_SeqAIJ, 1390 MatGetRow_SeqAIJ,MatRestoreRow_SeqAIJ, 1391 MatMult_SeqAIJ,MatMultAdd_SeqAIJ, 1392 MatMultTrans_SeqAIJ,MatMultTransAdd_SeqAIJ, 1393 MatSolve_SeqAIJ,MatSolveAdd_SeqAIJ, 1394 MatSolveTrans_SeqAIJ,MatSolveTransAdd_SeqAIJ, 1395 MatLUFactor_SeqAIJ,0, 1396 MatRelax_SeqAIJ, 1397 MatTranspose_SeqAIJ, 1398 MatGetInfo_SeqAIJ,MatEqual_SeqAIJ, 1399 MatGetDiagonal_SeqAIJ,MatDiagonalScale_SeqAIJ,MatNorm_SeqAIJ, 1400 0,MatAssemblyEnd_SeqAIJ, 1401 MatCompress_SeqAIJ, 1402 MatSetOption_SeqAIJ,MatZeroEntries_SeqAIJ,MatZeroRows_SeqAIJ, 1403 MatLUFactorSymbolic_SeqAIJ,MatLUFactorNumeric_SeqAIJ,0,0, 1404 MatGetSize_SeqAIJ,MatGetSize_SeqAIJ,MatGetOwnershipRange_SeqAIJ, 1405 MatILUFactorSymbolic_SeqAIJ,0, 1406 0,0,MatConvert_SeqAIJ, 1407 MatConvertSameType_SeqAIJ,0,0, 1408 MatILUFactor_SeqAIJ,0,0, 1409 MatGetSubMatrices_SeqAIJ,MatIncreaseOverlap_SeqAIJ, 1410 MatGetValues_SeqAIJ,0, 1411 MatPrintHelp_SeqAIJ, 1412 MatScale_SeqAIJ,0,0, 1413 MatILUDTFactor_SeqAIJ, 1414 MatGetBlockSize_SeqAIJ, 1415 MatGetRowIJ_SeqAIJ, 1416 MatRestoreRowIJ_SeqAIJ, 1417 MatGetColumnIJ_SeqAIJ, 1418 MatRestoreColumnIJ_SeqAIJ, 1419 MatFDColoringCreate_SeqAIJ, 1420 MatColoringPatch_SeqAIJ}; 1421 1422 extern int MatUseSuperLU_SeqAIJ(Mat); 1423 extern int MatUseEssl_SeqAIJ(Mat); 1424 extern int MatUseDXML_SeqAIJ(Mat); 1425 1426 /*@C 1427 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 1428 (the default parallel PETSc format). For good matrix assembly performance 1429 the user should preallocate the matrix storage by setting the parameter nz 1430 (or the array nzz). By setting these parameters accurately, performance 1431 during matrix assembly can be increased by more than a factor of 50. 1432 1433 Input Parameters: 1434 . comm - MPI communicator, set to MPI_COMM_SELF 1435 . m - number of rows 1436 . n - number of columns 1437 . nz - number of nonzeros per row (same for all rows) 1438 . nzz - array containing the number of nonzeros in the various rows 1439 (possibly different for each row) or PETSC_NULL 1440 1441 Output Parameter: 1442 . A - the matrix 1443 1444 Notes: 1445 The AIJ format (also called the Yale sparse matrix format or 1446 compressed row storage), is fully compatible with standard Fortran 77 1447 storage. That is, the stored row and column indices can begin at 1448 either one (as in Fortran) or zero. See the users' manual for details. 1449 1450 Specify the preallocated storage with either nz or nnz (not both). 1451 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1452 allocation. For large problems you MUST preallocate memory or you 1453 will get TERRIBLE performance, see the users' manual chapter on matrices. 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