1 2 3 #ifndef lint 4 static char vcid[] = "$Id: aij.c,v 1.198 1996/11/29 22:26:43 curfman Exp balay $"; 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->reallocs = 0; 548 A->info.nz_unneeded = (double)fshift; 549 550 /* check out for identical nodes. If found, use inode functions */ 551 ierr = Mat_AIJ_CheckInode(A); CHKERRQ(ierr); 552 return 0; 553 } 554 555 static int MatZeroEntries_SeqAIJ(Mat A) 556 { 557 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 558 PetscMemzero(a->a,(a->i[a->m]+a->indexshift)*sizeof(Scalar)); 559 return 0; 560 } 561 562 int MatDestroy_SeqAIJ(PetscObject obj) 563 { 564 Mat A = (Mat) obj; 565 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 566 int ierr; 567 568 #if defined(PETSC_LOG) 569 PLogObjectState(obj,"Rows=%d, Cols=%d, NZ=%d",a->m,a->n,a->nz); 570 #endif 571 PetscFree(a->a); 572 if (!a->singlemalloc) { PetscFree(a->i); PetscFree(a->j);} 573 if (a->diag) PetscFree(a->diag); 574 if (a->ilen) PetscFree(a->ilen); 575 if (a->imax) PetscFree(a->imax); 576 if (a->solve_work) PetscFree(a->solve_work); 577 if (a->inode.size) PetscFree(a->inode.size); 578 PetscFree(a); 579 if (A->mapping) { 580 ierr = ISLocalToGlobalMappingDestroy(A->mapping); CHKERRQ(ierr); 581 } 582 PLogObjectDestroy(A); 583 PetscHeaderDestroy(A); 584 return 0; 585 } 586 587 static int MatCompress_SeqAIJ(Mat A) 588 { 589 return 0; 590 } 591 592 static int MatSetOption_SeqAIJ(Mat A,MatOption op) 593 { 594 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 595 if (op == MAT_ROW_ORIENTED) a->roworiented = 1; 596 else if (op == MAT_COLUMN_ORIENTED) a->roworiented = 0; 597 else if (op == MAT_COLUMNS_SORTED) a->sorted = 1; 598 else if (op == MAT_COLUMNS_UNSORTED) a->sorted = 0; 599 else if (op == MAT_NO_NEW_NONZERO_LOCATIONS) a->nonew = 1; 600 else if (op == MAT_YES_NEW_NONZERO_LOCATIONS) a->nonew = 0; 601 else if (op == MAT_ROWS_SORTED || 602 op == MAT_ROWS_UNSORTED || 603 op == MAT_SYMMETRIC || 604 op == MAT_STRUCTURALLY_SYMMETRIC || 605 op == MAT_YES_NEW_DIAGONALS || 606 op == MAT_IGNORE_OFF_PROCESSOR_ENTRIES) 607 PLogInfo(A,"Info:MatSetOption_SeqAIJ:Option ignored\n"); 608 else if (op == MAT_NO_NEW_DIAGONALS) 609 {SETERRQ(PETSC_ERR_SUP,"MatSetOption_SeqAIJ:MAT_NO_NEW_DIAGONALS");} 610 else if (op == MAT_INODE_LIMIT_1) a->inode.limit = 1; 611 else if (op == MAT_INODE_LIMIT_2) a->inode.limit = 2; 612 else if (op == MAT_INODE_LIMIT_3) a->inode.limit = 3; 613 else if (op == MAT_INODE_LIMIT_4) a->inode.limit = 4; 614 else if (op == MAT_INODE_LIMIT_5) a->inode.limit = 5; 615 else 616 {SETERRQ(PETSC_ERR_SUP,"MatSetOption_SeqAIJ:unknown option");} 617 return 0; 618 } 619 620 static int MatGetDiagonal_SeqAIJ(Mat A,Vec v) 621 { 622 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 623 int i,j, n,shift = a->indexshift; 624 Scalar *x, zero = 0.0; 625 626 VecSet(&zero,v); 627 VecGetArray_Fast(v,x); VecGetLocalSize(v,&n); 628 if (n != a->m) SETERRQ(1,"MatGetDiagonal_SeqAIJ:Nonconforming matrix and vector"); 629 for ( i=0; i<a->m; i++ ) { 630 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 631 if (a->j[j]+shift == i) { 632 x[i] = a->a[j]; 633 break; 634 } 635 } 636 } 637 return 0; 638 } 639 640 /* -------------------------------------------------------*/ 641 /* Should check that shapes of vectors and matrices match */ 642 /* -------------------------------------------------------*/ 643 int MatMultTrans_SeqAIJ(Mat A,Vec xx,Vec yy) 644 { 645 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 646 Scalar *x, *y, *v, alpha; 647 int m = a->m, n, i, *idx, shift = a->indexshift; 648 649 VecGetArray_Fast(xx,x); VecGetArray_Fast(yy,y); 650 PetscMemzero(y,a->n*sizeof(Scalar)); 651 y = y + shift; /* shift for Fortran start by 1 indexing */ 652 for ( i=0; i<m; i++ ) { 653 idx = a->j + a->i[i] + shift; 654 v = a->a + a->i[i] + shift; 655 n = a->i[i+1] - a->i[i]; 656 alpha = x[i]; 657 while (n-->0) {y[*idx++] += alpha * *v++;} 658 } 659 PLogFlops(2*a->nz - a->n); 660 return 0; 661 } 662 663 int MatMultTransAdd_SeqAIJ(Mat A,Vec xx,Vec zz,Vec yy) 664 { 665 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 666 Scalar *x, *y, *v, alpha; 667 int m = a->m, n, i, *idx,shift = a->indexshift; 668 669 VecGetArray_Fast(xx,x); VecGetArray_Fast(yy,y); 670 if (zz != yy) VecCopy(zz,yy); 671 y = y + shift; /* shift for Fortran start by 1 indexing */ 672 for ( i=0; i<m; i++ ) { 673 idx = a->j + a->i[i] + shift; 674 v = a->a + a->i[i] + shift; 675 n = a->i[i+1] - a->i[i]; 676 alpha = x[i]; 677 while (n-->0) {y[*idx++] += alpha * *v++;} 678 } 679 PLogFlops(2*a->nz); 680 return 0; 681 } 682 683 int MatMult_SeqAIJ(Mat A,Vec xx,Vec yy) 684 { 685 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 686 Scalar *x, *y, *v, sum; 687 int m = a->m, n, i, *idx, shift = a->indexshift,*ii; 688 689 VecGetArray_Fast(xx,x); VecGetArray_Fast(yy,y); 690 x = x + shift; /* shift for Fortran start by 1 indexing */ 691 idx = a->j; 692 v = a->a; 693 ii = a->i; 694 for ( i=0; i<m; i++ ) { 695 n = ii[1] - ii[0]; ii++; 696 sum = 0.0; 697 /* SPARSEDENSEDOT(sum,x,v,idx,n); */ 698 /* for ( j=n-1; j>-1; j--) sum += v[j]*x[idx[j]]; */ 699 while (n--) sum += *v++ * x[*idx++]; 700 y[i] = sum; 701 } 702 PLogFlops(2*a->nz - m); 703 return 0; 704 } 705 706 int MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 707 { 708 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 709 Scalar *x, *y, *z, *v, sum; 710 int m = a->m, n, i, *idx, shift = a->indexshift,*ii; 711 712 VecGetArray_Fast(xx,x); VecGetArray_Fast(yy,y); VecGetArray_Fast(zz,z); 713 x = x + shift; /* shift for Fortran start by 1 indexing */ 714 idx = a->j; 715 v = a->a; 716 ii = a->i; 717 for ( i=0; i<m; i++ ) { 718 n = ii[1] - ii[0]; ii++; 719 sum = y[i]; 720 /* SPARSEDENSEDOT(sum,x,v,idx,n); */ 721 while (n--) sum += *v++ * x[*idx++]; 722 z[i] = sum; 723 } 724 PLogFlops(2*a->nz); 725 return 0; 726 } 727 728 /* 729 Adds diagonal pointers to sparse matrix structure. 730 */ 731 732 int MatMarkDiag_SeqAIJ(Mat A) 733 { 734 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 735 int i,j, *diag, m = a->m,shift = a->indexshift; 736 737 diag = (int *) PetscMalloc( (m+1)*sizeof(int)); CHKPTRQ(diag); 738 PLogObjectMemory(A,(m+1)*sizeof(int)); 739 for ( i=0; i<a->m; i++ ) { 740 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 741 if (a->j[j]+shift == i) { 742 diag[i] = j - shift; 743 break; 744 } 745 } 746 } 747 a->diag = diag; 748 return 0; 749 } 750 751 int MatRelax_SeqAIJ(Mat A,Vec bb,double omega,MatSORType flag, 752 double fshift,int its,Vec xx) 753 { 754 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 755 Scalar *x, *b, *bs, d, *xs, sum, *v = a->a,*t,scale,*ts, *xb; 756 int ierr, *idx, *diag,n = a->n, m = a->m, i, shift = a->indexshift; 757 758 VecGetArray_Fast(xx,x); VecGetArray_Fast(bb,b); 759 if (!a->diag) {if ((ierr = MatMarkDiag_SeqAIJ(A))) return ierr;} 760 diag = a->diag; 761 xs = x + shift; /* shifted by one for index start of a or a->j*/ 762 if (flag == SOR_APPLY_UPPER) { 763 /* apply ( U + D/omega) to the vector */ 764 bs = b + shift; 765 for ( i=0; i<m; i++ ) { 766 d = fshift + a->a[diag[i] + shift]; 767 n = a->i[i+1] - diag[i] - 1; 768 idx = a->j + diag[i] + (!shift); 769 v = a->a + diag[i] + (!shift); 770 sum = b[i]*d/omega; 771 SPARSEDENSEDOT(sum,bs,v,idx,n); 772 x[i] = sum; 773 } 774 return 0; 775 } 776 if (flag == SOR_APPLY_LOWER) { 777 SETERRQ(1,"MatRelax_SeqAIJ:SOR_APPLY_LOWER is not done"); 778 } 779 else if (flag & SOR_EISENSTAT) { 780 /* Let A = L + U + D; where L is lower trianglar, 781 U is upper triangular, E is diagonal; This routine applies 782 783 (L + E)^{-1} A (U + E)^{-1} 784 785 to a vector efficiently using Eisenstat's trick. This is for 786 the case of SSOR preconditioner, so E is D/omega where omega 787 is the relaxation factor. 788 */ 789 t = (Scalar *) PetscMalloc( m*sizeof(Scalar) ); CHKPTRQ(t); 790 scale = (2.0/omega) - 1.0; 791 792 /* x = (E + U)^{-1} b */ 793 for ( i=m-1; i>=0; i-- ) { 794 d = fshift + a->a[diag[i] + shift]; 795 n = a->i[i+1] - diag[i] - 1; 796 idx = a->j + diag[i] + (!shift); 797 v = a->a + diag[i] + (!shift); 798 sum = b[i]; 799 SPARSEDENSEMDOT(sum,xs,v,idx,n); 800 x[i] = omega*(sum/d); 801 } 802 803 /* t = b - (2*E - D)x */ 804 v = a->a; 805 for ( i=0; i<m; i++ ) { t[i] = b[i] - scale*(v[*diag++ + shift])*x[i]; } 806 807 /* t = (E + L)^{-1}t */ 808 ts = t + shift; /* shifted by one for index start of a or a->j*/ 809 diag = a->diag; 810 for ( i=0; i<m; i++ ) { 811 d = fshift + a->a[diag[i]+shift]; 812 n = diag[i] - a->i[i]; 813 idx = a->j + a->i[i] + shift; 814 v = a->a + a->i[i] + shift; 815 sum = t[i]; 816 SPARSEDENSEMDOT(sum,ts,v,idx,n); 817 t[i] = omega*(sum/d); 818 } 819 820 /* x = x + t */ 821 for ( i=0; i<m; i++ ) { x[i] += t[i]; } 822 PetscFree(t); 823 return 0; 824 } 825 if (flag & SOR_ZERO_INITIAL_GUESS) { 826 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 827 for ( i=0; i<m; i++ ) { 828 d = fshift + a->a[diag[i]+shift]; 829 n = diag[i] - a->i[i]; 830 idx = a->j + a->i[i] + shift; 831 v = a->a + a->i[i] + shift; 832 sum = b[i]; 833 SPARSEDENSEMDOT(sum,xs,v,idx,n); 834 x[i] = omega*(sum/d); 835 } 836 xb = x; 837 } 838 else xb = b; 839 if ((flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) && 840 (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP)) { 841 for ( i=0; i<m; i++ ) { 842 x[i] *= a->a[diag[i]+shift]; 843 } 844 } 845 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 846 for ( i=m-1; i>=0; i-- ) { 847 d = fshift + a->a[diag[i] + shift]; 848 n = a->i[i+1] - diag[i] - 1; 849 idx = a->j + diag[i] + (!shift); 850 v = a->a + diag[i] + (!shift); 851 sum = xb[i]; 852 SPARSEDENSEMDOT(sum,xs,v,idx,n); 853 x[i] = omega*(sum/d); 854 } 855 } 856 its--; 857 } 858 while (its--) { 859 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 860 for ( i=0; i<m; i++ ) { 861 d = fshift + a->a[diag[i]+shift]; 862 n = a->i[i+1] - a->i[i]; 863 idx = a->j + a->i[i] + shift; 864 v = a->a + a->i[i] + shift; 865 sum = b[i]; 866 SPARSEDENSEMDOT(sum,xs,v,idx,n); 867 x[i] = (1. - omega)*x[i] + omega*(sum + a->a[diag[i]+shift]*x[i])/d; 868 } 869 } 870 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 871 for ( i=m-1; i>=0; i-- ) { 872 d = fshift + a->a[diag[i] + shift]; 873 n = a->i[i+1] - a->i[i]; 874 idx = a->j + a->i[i] + shift; 875 v = a->a + a->i[i] + shift; 876 sum = b[i]; 877 SPARSEDENSEMDOT(sum,xs,v,idx,n); 878 x[i] = (1. - omega)*x[i] + omega*(sum + a->a[diag[i]+shift]*x[i])/d; 879 } 880 } 881 } 882 return 0; 883 } 884 885 static int MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 886 { 887 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 888 889 info->rows_global = (double)a->m; 890 info->columns_global = (double)a->n; 891 info->rows_local = (double)a->m; 892 info->columns_local = (double)a->n; 893 info->block_size = 1.0; 894 info->nz_allocated = (double)a->maxnz; 895 info->nz_used = (double)a->nz; 896 info->nz_unneeded = (double)(a->maxnz - a->nz); 897 /* if (info->nz_unneeded != A->info.nz_unneeded) 898 printf("space descrepancy: maxnz-nz = %d, nz_unneeded = %d\n",(int)info->nz_unneeded,(int)A->info.nz_unneeded); */ 899 info->assemblies = (double)A->num_ass; 900 info->mallocs = (double)a->reallocs; 901 info->memory = A->mem; 902 if (A->factor) { 903 info->fill_ratio_given = A->info.fill_ratio_given; 904 info->fill_ratio_needed = A->info.fill_ratio_needed; 905 info->factor_mallocs = A->info.factor_mallocs; 906 } else { 907 info->fill_ratio_given = 0; 908 info->fill_ratio_needed = 0; 909 info->factor_mallocs = 0; 910 } 911 return 0; 912 } 913 914 extern int MatLUFactorSymbolic_SeqAIJ(Mat,IS,IS,double,Mat*); 915 extern int MatLUFactorNumeric_SeqAIJ(Mat,Mat*); 916 extern int MatLUFactor_SeqAIJ(Mat,IS,IS,double); 917 extern int MatSolve_SeqAIJ(Mat,Vec,Vec); 918 extern int MatSolveAdd_SeqAIJ(Mat,Vec,Vec,Vec); 919 extern int MatSolveTrans_SeqAIJ(Mat,Vec,Vec); 920 extern int MatSolveTransAdd_SeqAIJ(Mat,Vec,Vec,Vec); 921 922 static int MatZeroRows_SeqAIJ(Mat A,IS is,Scalar *diag) 923 { 924 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 925 int i,ierr,N, *rows,m = a->m - 1,shift = a->indexshift; 926 927 ierr = ISGetSize(is,&N); CHKERRQ(ierr); 928 ierr = ISGetIndices(is,&rows); CHKERRQ(ierr); 929 if (diag) { 930 for ( i=0; i<N; i++ ) { 931 if (rows[i] < 0 || rows[i] > m) SETERRQ(1,"MatZeroRows_SeqAIJ:row out of range"); 932 if (a->ilen[rows[i]] > 0) { /* in case row was completely empty */ 933 a->ilen[rows[i]] = 1; 934 a->a[a->i[rows[i]]+shift] = *diag; 935 a->j[a->i[rows[i]]+shift] = rows[i]+shift; 936 } 937 else { 938 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],diag,INSERT_VALUES); 939 CHKERRQ(ierr); 940 } 941 } 942 } 943 else { 944 for ( i=0; i<N; i++ ) { 945 if (rows[i] < 0 || rows[i] > m) SETERRQ(1,"MatZeroRows_SeqAIJ:row out of range"); 946 a->ilen[rows[i]] = 0; 947 } 948 } 949 ISRestoreIndices(is,&rows); 950 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 951 return 0; 952 } 953 954 static int MatGetSize_SeqAIJ(Mat A,int *m,int *n) 955 { 956 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 957 *m = a->m; *n = a->n; 958 return 0; 959 } 960 961 static int MatGetOwnershipRange_SeqAIJ(Mat A,int *m,int *n) 962 { 963 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 964 *m = 0; *n = a->m; 965 return 0; 966 } 967 int MatGetRow_SeqAIJ(Mat A,int row,int *nz,int **idx,Scalar **v) 968 { 969 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 970 int *itmp,i,shift = a->indexshift; 971 972 if (row < 0 || row >= a->m) SETERRQ(1,"MatGetRow_SeqAIJ:Row out of range"); 973 974 *nz = a->i[row+1] - a->i[row]; 975 if (v) *v = a->a + a->i[row] + shift; 976 if (idx) { 977 itmp = a->j + a->i[row] + shift; 978 if (*nz && shift) { 979 *idx = (int *) PetscMalloc( (*nz)*sizeof(int) ); CHKPTRQ(*idx); 980 for ( i=0; i<(*nz); i++ ) {(*idx)[i] = itmp[i] + shift;} 981 } else if (*nz) { 982 *idx = itmp; 983 } 984 else *idx = 0; 985 } 986 return 0; 987 } 988 989 int MatRestoreRow_SeqAIJ(Mat A,int row,int *nz,int **idx,Scalar **v) 990 { 991 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 992 if (idx) {if (*idx && a->indexshift) PetscFree(*idx);} 993 return 0; 994 } 995 996 static int MatNorm_SeqAIJ(Mat A,NormType type,double *norm) 997 { 998 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 999 Scalar *v = a->a; 1000 double sum = 0.0; 1001 int i, j,shift = a->indexshift; 1002 1003 if (type == NORM_FROBENIUS) { 1004 for (i=0; i<a->nz; i++ ) { 1005 #if defined(PETSC_COMPLEX) 1006 sum += real(conj(*v)*(*v)); v++; 1007 #else 1008 sum += (*v)*(*v); v++; 1009 #endif 1010 } 1011 *norm = sqrt(sum); 1012 } 1013 else if (type == NORM_1) { 1014 double *tmp; 1015 int *jj = a->j; 1016 tmp = (double *) PetscMalloc( a->n*sizeof(double) ); CHKPTRQ(tmp); 1017 PetscMemzero(tmp,a->n*sizeof(double)); 1018 *norm = 0.0; 1019 for ( j=0; j<a->nz; j++ ) { 1020 tmp[*jj++ + shift] += PetscAbsScalar(*v); v++; 1021 } 1022 for ( j=0; j<a->n; j++ ) { 1023 if (tmp[j] > *norm) *norm = tmp[j]; 1024 } 1025 PetscFree(tmp); 1026 } 1027 else if (type == NORM_INFINITY) { 1028 *norm = 0.0; 1029 for ( j=0; j<a->m; j++ ) { 1030 v = a->a + a->i[j] + shift; 1031 sum = 0.0; 1032 for ( i=0; i<a->i[j+1]-a->i[j]; i++ ) { 1033 sum += PetscAbsScalar(*v); v++; 1034 } 1035 if (sum > *norm) *norm = sum; 1036 } 1037 } 1038 else { 1039 SETERRQ(1,"MatNorm_SeqAIJ:No support for two norm yet"); 1040 } 1041 return 0; 1042 } 1043 1044 static int MatTranspose_SeqAIJ(Mat A,Mat *B) 1045 { 1046 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1047 Mat C; 1048 int i, ierr, *aj = a->j, *ai = a->i, m = a->m, len, *col; 1049 int shift = a->indexshift; 1050 Scalar *array = a->a; 1051 1052 if (B == PETSC_NULL && m != a->n) 1053 SETERRQ(1,"MatTranspose_SeqAIJ:Square matrix only for in-place"); 1054 col = (int *) PetscMalloc((1+a->n)*sizeof(int)); CHKPTRQ(col); 1055 PetscMemzero(col,(1+a->n)*sizeof(int)); 1056 if (shift) { 1057 for ( i=0; i<ai[m]-1; i++ ) aj[i] -= 1; 1058 } 1059 for ( i=0; i<ai[m]+shift; i++ ) col[aj[i]] += 1; 1060 ierr = MatCreateSeqAIJ(A->comm,a->n,m,0,col,&C); CHKERRQ(ierr); 1061 PetscFree(col); 1062 for ( i=0; i<m; i++ ) { 1063 len = ai[i+1]-ai[i]; 1064 ierr = MatSetValues(C,len,aj,1,&i,array,INSERT_VALUES); CHKERRQ(ierr); 1065 array += len; aj += len; 1066 } 1067 if (shift) { 1068 for ( i=0; i<ai[m]-1; i++ ) aj[i] += 1; 1069 } 1070 1071 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1072 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1073 1074 if (B != PETSC_NULL) { 1075 *B = C; 1076 } else { 1077 /* This isn't really an in-place transpose */ 1078 PetscFree(a->a); 1079 if (!a->singlemalloc) {PetscFree(a->i); PetscFree(a->j);} 1080 if (a->diag) PetscFree(a->diag); 1081 if (a->ilen) PetscFree(a->ilen); 1082 if (a->imax) PetscFree(a->imax); 1083 if (a->solve_work) PetscFree(a->solve_work); 1084 if (a->inode.size) PetscFree(a->inode.size); 1085 PetscFree(a); 1086 PetscMemcpy(A,C,sizeof(struct _Mat)); 1087 PetscHeaderDestroy(C); 1088 } 1089 return 0; 1090 } 1091 1092 static int MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 1093 { 1094 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1095 Scalar *l,*r,x,*v; 1096 int i,j,m = a->m, n = a->n, M, nz = a->nz, *jj,shift = a->indexshift; 1097 1098 if (ll) { 1099 /* The local size is used so that VecMPI can be passed to this routine 1100 by MatDiagonalScale_MPIAIJ */ 1101 VecGetLocalSize_Fast(ll,m); 1102 if (m != a->m) SETERRQ(1,"MatDiagonalScale_SeqAIJ:Left scaling vector wrong length"); 1103 VecGetArray_Fast(ll,l); 1104 v = a->a; 1105 for ( i=0; i<m; i++ ) { 1106 x = l[i]; 1107 M = a->i[i+1] - a->i[i]; 1108 for ( j=0; j<M; j++ ) { (*v++) *= x;} 1109 } 1110 PLogFlops(nz); 1111 } 1112 if (rr) { 1113 VecGetLocalSize_Fast(rr,n); 1114 if (n != a->n) SETERRQ(1,"MatDiagonalScale_SeqAIJ:Right scaling vector wrong length"); 1115 VecGetArray_Fast(rr,r); 1116 v = a->a; jj = a->j; 1117 for ( i=0; i<nz; i++ ) { 1118 (*v++) *= r[*jj++ + shift]; 1119 } 1120 PLogFlops(nz); 1121 } 1122 return 0; 1123 } 1124 1125 static int MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,MatGetSubMatrixCall scall,Mat *B) 1126 { 1127 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data,*c; 1128 int nznew, *smap, i, k, kstart, kend, ierr, oldcols = a->n,*lens; 1129 int row,mat_i,*mat_j,tcol,first,step,*mat_ilen; 1130 register int sum,lensi; 1131 int *irow, *icol, nrows, ncols, shift = a->indexshift,*ssmap; 1132 int *starts,*j_new,*i_new,*aj = a->j, *ai = a->i,ii,*ailen = a->ilen; 1133 Scalar *a_new,*mat_a; 1134 Mat C; 1135 1136 ierr = ISSorted(isrow,(PetscTruth*)&i); 1137 if (!i) SETERRQ(1,"MatGetSubmatrices_SeqAIJ:ISrow is not sorted"); 1138 ierr = ISSorted(iscol,(PetscTruth*)&i); 1139 if (!i) SETERRQ(1,"MatGetSubmatrices_SeqAIJ:IScol is not sorted"); 1140 1141 ierr = ISGetIndices(isrow,&irow); CHKERRQ(ierr); 1142 ierr = ISGetSize(isrow,&nrows); CHKERRQ(ierr); 1143 ierr = ISGetSize(iscol,&ncols); CHKERRQ(ierr); 1144 1145 if (ISStrideGetInfo(iscol,&first,&step) && step == 1) { /* no need to sort */ 1146 /* special case of contiguous rows */ 1147 lens = (int *) PetscMalloc((ncols+nrows+1)*sizeof(int)); CHKPTRQ(lens); 1148 starts = lens + ncols; 1149 /* loop over new rows determining lens and starting points */ 1150 for (i=0; i<nrows; i++) { 1151 kstart = ai[irow[i]]+shift; 1152 kend = kstart + ailen[irow[i]]; 1153 for ( k=kstart; k<kend; k++ ) { 1154 if (aj[k]+shift >= first) { 1155 starts[i] = k; 1156 break; 1157 } 1158 } 1159 sum = 0; 1160 while (k < kend) { 1161 if (aj[k++]+shift >= first+ncols) break; 1162 sum++; 1163 } 1164 lens[i] = sum; 1165 } 1166 /* create submatrix */ 1167 if (scall == MAT_REUSE_MATRIX) { 1168 int n_cols,n_rows; 1169 ierr = MatGetSize(*B,&n_rows,&n_cols); CHKERRQ(ierr); 1170 if (n_rows != nrows || n_cols != ncols) SETERRQ(1,"MatGetSubMatrix_SeqAIJ:"); 1171 ierr = MatZeroEntries(*B); CHKERRQ(ierr); 1172 C = *B; 1173 } 1174 else { 1175 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1176 } 1177 c = (Mat_SeqAIJ*) C->data; 1178 1179 /* loop over rows inserting into submatrix */ 1180 a_new = c->a; 1181 j_new = c->j; 1182 i_new = c->i; 1183 i_new[0] = -shift; 1184 for (i=0; i<nrows; i++) { 1185 ii = starts[i]; 1186 lensi = lens[i]; 1187 for ( k=0; k<lensi; k++ ) { 1188 *j_new++ = aj[ii+k] - first; 1189 } 1190 PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(Scalar)); 1191 a_new += lensi; 1192 i_new[i+1] = i_new[i] + lensi; 1193 c->ilen[i] = lensi; 1194 } 1195 PetscFree(lens); 1196 } 1197 else { 1198 ierr = ISGetIndices(iscol,&icol); CHKERRQ(ierr); 1199 smap = (int *) PetscMalloc((1+oldcols)*sizeof(int)); CHKPTRQ(smap); 1200 ssmap = smap + shift; 1201 lens = (int *) PetscMalloc((1+nrows)*sizeof(int)); CHKPTRQ(lens); 1202 PetscMemzero(smap,oldcols*sizeof(int)); 1203 for ( i=0; i<ncols; i++ ) smap[icol[i]] = i+1; 1204 /* determine lens of each row */ 1205 for (i=0; i<nrows; i++) { 1206 kstart = ai[irow[i]]+shift; 1207 kend = kstart + a->ilen[irow[i]]; 1208 lens[i] = 0; 1209 for ( k=kstart; k<kend; k++ ) { 1210 if (ssmap[aj[k]]) { 1211 lens[i]++; 1212 } 1213 } 1214 } 1215 /* Create and fill new matrix */ 1216 if (scall == MAT_REUSE_MATRIX) { 1217 c = (Mat_SeqAIJ *)((*B)->data); 1218 1219 if (c->m != nrows || c->n != ncols) SETERRQ(1,"MatGetSubMatrix_SeqAIJ:"); 1220 if (PetscMemcmp(c->ilen,lens, c->m *sizeof(int))) { 1221 SETERRQ(1,"MatGetSubmatrices_SeqAIJ:Cannot reuse matrix. wrong no of nonzeros"); 1222 } 1223 PetscMemzero(c->ilen,c->m*sizeof(int)); 1224 C = *B; 1225 } 1226 else { 1227 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1228 } 1229 c = (Mat_SeqAIJ *)(C->data); 1230 for (i=0; i<nrows; i++) { 1231 row = irow[i]; 1232 nznew = 0; 1233 kstart = ai[row]+shift; 1234 kend = kstart + a->ilen[row]; 1235 mat_i = c->i[i]+shift; 1236 mat_j = c->j + mat_i; 1237 mat_a = c->a + mat_i; 1238 mat_ilen = c->ilen + i; 1239 for ( k=kstart; k<kend; k++ ) { 1240 if ((tcol=ssmap[a->j[k]])) { 1241 *mat_j++ = tcol - (!shift); 1242 *mat_a++ = a->a[k]; 1243 (*mat_ilen)++; 1244 1245 } 1246 } 1247 } 1248 /* Free work space */ 1249 ierr = ISRestoreIndices(iscol,&icol); CHKERRQ(ierr); 1250 PetscFree(smap); PetscFree(lens); 1251 } 1252 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1253 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1254 1255 ierr = ISRestoreIndices(isrow,&irow); CHKERRQ(ierr); 1256 *B = C; 1257 return 0; 1258 } 1259 1260 /* 1261 note: This can only work for identity for row and col. It would 1262 be good to check this and otherwise generate an error. 1263 */ 1264 static int MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,double efill,int fill) 1265 { 1266 Mat_SeqAIJ *a = (Mat_SeqAIJ *) inA->data; 1267 int ierr; 1268 Mat outA; 1269 1270 if (fill != 0) SETERRQ(1,"MatILUFactor_SeqAIJ:Only fill=0 supported"); 1271 1272 outA = inA; 1273 inA->factor = FACTOR_LU; 1274 a->row = row; 1275 a->col = col; 1276 1277 a->solve_work = (Scalar *) PetscMalloc( (a->m+1)*sizeof(Scalar)); CHKPTRQ(a->solve_work); 1278 1279 if (!a->diag) { 1280 ierr = MatMarkDiag_SeqAIJ(inA); CHKERRQ(ierr); 1281 } 1282 ierr = MatLUFactorNumeric_SeqAIJ(inA,&outA); CHKERRQ(ierr); 1283 return 0; 1284 } 1285 1286 #include "pinclude/plapack.h" 1287 static int MatScale_SeqAIJ(Scalar *alpha,Mat inA) 1288 { 1289 Mat_SeqAIJ *a = (Mat_SeqAIJ *) inA->data; 1290 int one = 1; 1291 BLscal_( &a->nz, alpha, a->a, &one ); 1292 PLogFlops(a->nz); 1293 return 0; 1294 } 1295 1296 static int MatGetSubMatrices_SeqAIJ(Mat A,int n, IS *irow,IS *icol,MatGetSubMatrixCall scall, 1297 Mat **B) 1298 { 1299 int ierr,i; 1300 1301 if (scall == MAT_INITIAL_MATRIX) { 1302 *B = (Mat *) PetscMalloc( (n+1)*sizeof(Mat) ); CHKPTRQ(*B); 1303 } 1304 1305 for ( i=0; i<n; i++ ) { 1306 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],scall,&(*B)[i]);CHKERRQ(ierr); 1307 } 1308 return 0; 1309 } 1310 1311 static int MatGetBlockSize_SeqAIJ(Mat A, int *bs) 1312 { 1313 *bs = 1; 1314 return 0; 1315 } 1316 1317 static int MatIncreaseOverlap_SeqAIJ(Mat A, int is_max, IS *is, int ov) 1318 { 1319 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1320 int shift, row, i,j,k,l,m,n, *idx,ierr, *nidx, isz, val; 1321 int start, end, *ai, *aj; 1322 char *table; 1323 shift = a->indexshift; 1324 m = a->m; 1325 ai = a->i; 1326 aj = a->j+shift; 1327 1328 if (ov < 0) SETERRQ(1,"MatIncreaseOverlap_SeqAIJ: illegal overlap value used"); 1329 1330 table = (char *) PetscMalloc((m/BITSPERBYTE +1)*sizeof(char)); CHKPTRQ(table); 1331 nidx = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(nidx); 1332 1333 for ( i=0; i<is_max; i++ ) { 1334 /* Initialise the two local arrays */ 1335 isz = 0; 1336 PetscMemzero(table,(m/BITSPERBYTE +1)*sizeof(char)); 1337 1338 /* Extract the indices, assume there can be duplicate entries */ 1339 ierr = ISGetIndices(is[i],&idx); CHKERRQ(ierr); 1340 ierr = ISGetSize(is[i],&n); CHKERRQ(ierr); 1341 1342 /* Enter these into the temp arrays i.e mark table[row], enter row into new index */ 1343 for ( j=0; j<n ; ++j){ 1344 if(!BT_LOOKUP(table, idx[j])) { nidx[isz++] = idx[j];} 1345 } 1346 ierr = ISRestoreIndices(is[i],&idx); CHKERRQ(ierr); 1347 ierr = ISDestroy(is[i]); CHKERRQ(ierr); 1348 1349 k = 0; 1350 for ( j=0; j<ov; j++){ /* for each overlap*/ 1351 n = isz; 1352 for ( ; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 1353 row = nidx[k]; 1354 start = ai[row]; 1355 end = ai[row+1]; 1356 for ( l = start; l<end ; l++){ 1357 val = aj[l] + shift; 1358 if (!BT_LOOKUP(table,val)) {nidx[isz++] = val;} 1359 } 1360 } 1361 } 1362 ierr = ISCreateGeneral(MPI_COMM_SELF, isz, nidx, (is+i)); CHKERRQ(ierr); 1363 } 1364 PetscFree(table); 1365 PetscFree(nidx); 1366 return 0; 1367 } 1368 1369 int MatPrintHelp_SeqAIJ(Mat A) 1370 { 1371 static int called = 0; 1372 MPI_Comm comm = A->comm; 1373 1374 if (called) return 0; else called = 1; 1375 PetscPrintf(comm," Options for MATSEQAIJ and MATMPIAIJ matrix formats (the defaults):\n"); 1376 PetscPrintf(comm," -mat_lu_pivotthreshold <threshold>\n"); 1377 PetscPrintf(comm," -mat_aij_oneindex - internal indices begin at 1 instead of the default 0.\n"); 1378 PetscPrintf(comm," -mat_aij_no_inode - Do not use inodes\n"); 1379 PetscPrintf(comm," -mat_aij_inode_limit <limit> - Set inode limit (max limit=5)\n"); 1380 #if defined(HAVE_ESSL) 1381 PetscPrintf(comm," -mat_aij_essl - Use IBM sparse LU factorization and solve.\n"); 1382 #endif 1383 return 0; 1384 } 1385 static int MatEqual_SeqAIJ(Mat A,Mat B, PetscTruth* flg); 1386 extern int MatFDColoringCreate_SeqAIJ(Mat,ISColoring,MatFDColoring); 1387 extern int MatColoringPatch_SeqAIJ(Mat,int,int *,ISColoring *); 1388 1389 /* -------------------------------------------------------------------*/ 1390 static struct _MatOps MatOps = {MatSetValues_SeqAIJ, 1391 MatGetRow_SeqAIJ,MatRestoreRow_SeqAIJ, 1392 MatMult_SeqAIJ,MatMultAdd_SeqAIJ, 1393 MatMultTrans_SeqAIJ,MatMultTransAdd_SeqAIJ, 1394 MatSolve_SeqAIJ,MatSolveAdd_SeqAIJ, 1395 MatSolveTrans_SeqAIJ,MatSolveTransAdd_SeqAIJ, 1396 MatLUFactor_SeqAIJ,0, 1397 MatRelax_SeqAIJ, 1398 MatTranspose_SeqAIJ, 1399 MatGetInfo_SeqAIJ,MatEqual_SeqAIJ, 1400 MatGetDiagonal_SeqAIJ,MatDiagonalScale_SeqAIJ,MatNorm_SeqAIJ, 1401 0,MatAssemblyEnd_SeqAIJ, 1402 MatCompress_SeqAIJ, 1403 MatSetOption_SeqAIJ,MatZeroEntries_SeqAIJ,MatZeroRows_SeqAIJ, 1404 MatLUFactorSymbolic_SeqAIJ,MatLUFactorNumeric_SeqAIJ,0,0, 1405 MatGetSize_SeqAIJ,MatGetSize_SeqAIJ,MatGetOwnershipRange_SeqAIJ, 1406 MatILUFactorSymbolic_SeqAIJ,0, 1407 0,0,MatConvert_SeqAIJ, 1408 MatConvertSameType_SeqAIJ,0,0, 1409 MatILUFactor_SeqAIJ,0,0, 1410 MatGetSubMatrices_SeqAIJ,MatIncreaseOverlap_SeqAIJ, 1411 MatGetValues_SeqAIJ,0, 1412 MatPrintHelp_SeqAIJ, 1413 MatScale_SeqAIJ,0,0, 1414 MatILUDTFactor_SeqAIJ, 1415 MatGetBlockSize_SeqAIJ, 1416 MatGetRowIJ_SeqAIJ, 1417 MatRestoreRowIJ_SeqAIJ, 1418 MatGetColumnIJ_SeqAIJ, 1419 MatRestoreColumnIJ_SeqAIJ, 1420 MatFDColoringCreate_SeqAIJ, 1421 MatColoringPatch_SeqAIJ}; 1422 1423 extern int MatUseSuperLU_SeqAIJ(Mat); 1424 extern int MatUseEssl_SeqAIJ(Mat); 1425 extern int MatUseDXML_SeqAIJ(Mat); 1426 1427 /*@C 1428 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 1429 (the default parallel PETSc format). For good matrix assembly performance 1430 the user should preallocate the matrix storage by setting the parameter nz 1431 (or the array nzz). By setting these parameters accurately, performance 1432 during matrix assembly can be increased by more than a factor of 50. 1433 1434 Input Parameters: 1435 . comm - MPI communicator, set to MPI_COMM_SELF 1436 . m - number of rows 1437 . n - number of columns 1438 . nz - number of nonzeros per row (same for all rows) 1439 . nzz - array containing the number of nonzeros in the various rows 1440 (possibly different for each row) or PETSC_NULL 1441 1442 Output Parameter: 1443 . A - the matrix 1444 1445 Notes: 1446 The AIJ format (also called the Yale sparse matrix format or 1447 compressed row storage), is fully compatible with standard Fortran 77 1448 storage. That is, the stored row and column indices can begin at 1449 either one (as in Fortran) or zero. See the users' manual for details. 1450 1451 Specify the preallocated storage with either nz or nnz (not both). 1452 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1453 allocation. For large problems you MUST preallocate memory or you 1454 will get TERRIBLE performance, see the users' manual chapter on matrices. 1455 1456 By default, this format uses inodes (identical nodes) when possible, to 1457 improve numerical efficiency of Matrix vector products and solves. We 1458 search for consecutive rows with the same nonzero structure, thereby 1459 reusing matrix information to achieve increased efficiency. 1460 1461 Options Database Keys: 1462 $ -mat_aij_no_inode - Do not use inodes 1463 $ -mat_aij_inode_limit <limit> - Set inode limit. 1464 $ (max limit=5) 1465 $ -mat_aij_oneindex - Internally use indexing starting at 1 1466 $ rather than 0. Note: When calling MatSetValues(), 1467 $ the user still MUST index entries starting at 0! 1468 1469 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues() 1470 @*/ 1471 int MatCreateSeqAIJ(MPI_Comm comm,int m,int n,int nz,int *nnz, Mat *A) 1472 { 1473 Mat B; 1474 Mat_SeqAIJ *b; 1475 int i, len, ierr, flg,size; 1476 1477 MPI_Comm_size(comm,&size); 1478 if (size > 1) SETERRQ(1,"MatCreateSeqAIJ:Comm must be of size 1"); 1479 1480 *A = 0; 1481 PetscHeaderCreate(B,_Mat,MAT_COOKIE,MATSEQAIJ,comm); 1482 PLogObjectCreate(B); 1483 B->data = (void *) (b = PetscNew(Mat_SeqAIJ)); CHKPTRQ(b); 1484 PetscMemzero(b,sizeof(Mat_SeqAIJ)); 1485 PetscMemcpy(&B->ops,&MatOps,sizeof(struct _MatOps)); 1486 B->destroy = MatDestroy_SeqAIJ; 1487 B->view = MatView_SeqAIJ; 1488 B->factor = 0; 1489 B->lupivotthreshold = 1.0; 1490 B->mapping = 0; 1491 ierr = OptionsGetDouble(PETSC_NULL,"-mat_lu_pivotthreshold",&B->lupivotthreshold, 1492 &flg); CHKERRQ(ierr); 1493 b->ilu_preserve_row_sums = PETSC_FALSE; 1494 ierr = OptionsHasName(PETSC_NULL,"-pc_ilu_preserve_row_sums", 1495 (int*) &b->ilu_preserve_row_sums); CHKERRQ(ierr); 1496 b->row = 0; 1497 b->col = 0; 1498 b->indexshift = 0; 1499 b->reallocs = 0; 1500 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_oneindex", &flg); CHKERRQ(ierr); 1501 if (flg) b->indexshift = -1; 1502 1503 b->m = m; B->m = m; B->M = m; 1504 b->n = n; B->n = n; B->N = n; 1505 b->imax = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(b->imax); 1506 if (nnz == PETSC_NULL) { 1507 if (nz == PETSC_DEFAULT) nz = 10; 1508 else if (nz <= 0) nz = 1; 1509 for ( i=0; i<m; i++ ) b->imax[i] = nz; 1510 nz = nz*m; 1511 } 1512 else { 1513 nz = 0; 1514 for ( i=0; i<m; i++ ) {b->imax[i] = nnz[i]; nz += nnz[i];} 1515 } 1516 1517 /* allocate the matrix space */ 1518 len = nz*(sizeof(int) + sizeof(Scalar)) + (b->m+1)*sizeof(int); 1519 b->a = (Scalar *) PetscMalloc( len ); CHKPTRQ(b->a); 1520 b->j = (int *) (b->a + nz); 1521 PetscMemzero(b->j,nz*sizeof(int)); 1522 b->i = b->j + nz; 1523 b->singlemalloc = 1; 1524 1525 b->i[0] = -b->indexshift; 1526 for (i=1; i<m+1; i++) { 1527 b->i[i] = b->i[i-1] + b->imax[i-1]; 1528 } 1529 1530 /* b->ilen will count nonzeros in each row so far. */ 1531 b->ilen = (int *) PetscMalloc((m+1)*sizeof(int)); 1532 PLogObjectMemory(B,len+2*(m+1)*sizeof(int)+sizeof(struct _Mat)+sizeof(Mat_SeqAIJ)); 1533 for ( i=0; i<b->m; i++ ) { b->ilen[i] = 0;} 1534 1535 b->nz = 0; 1536 b->maxnz = nz; 1537 b->sorted = 0; 1538 b->roworiented = 1; 1539 b->nonew = 0; 1540 b->diag = 0; 1541 b->solve_work = 0; 1542 b->spptr = 0; 1543 b->inode.node_count = 0; 1544 b->inode.size = 0; 1545 b->inode.limit = 5; 1546 b->inode.max_limit = 5; 1547 B->info.nz_unneeded = (double)b->maxnz; 1548 1549 *A = B; 1550 1551 /* SuperLU is not currently supported through PETSc */ 1552 #if defined(HAVE_SUPERLU) 1553 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_superlu", &flg); CHKERRQ(ierr); 1554 if (flg) { ierr = MatUseSuperLU_SeqAIJ(B); CHKERRQ(ierr); } 1555 #endif 1556 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_essl", &flg); CHKERRQ(ierr); 1557 if (flg) { ierr = MatUseEssl_SeqAIJ(B); CHKERRQ(ierr); } 1558 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_dxml", &flg); CHKERRQ(ierr); 1559 if (flg) { 1560 if (!b->indexshift) SETERRQ(1,"MatCreateSeqAIJ:need -mat_aij_oneindex with -mat_aij_dxml"); 1561 ierr = MatUseDXML_SeqAIJ(B); CHKERRQ(ierr); 1562 } 1563 ierr = OptionsHasName(PETSC_NULL,"-help", &flg); CHKERRQ(ierr); 1564 if (flg) {ierr = MatPrintHelp(B); CHKERRQ(ierr); } 1565 return 0; 1566 } 1567 1568 int MatConvertSameType_SeqAIJ(Mat A,Mat *B,int cpvalues) 1569 { 1570 Mat C; 1571 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ *) A->data; 1572 int i,len, m = a->m,shift = a->indexshift; 1573 1574 *B = 0; 1575 PetscHeaderCreate(C,_Mat,MAT_COOKIE,MATSEQAIJ,A->comm); 1576 PLogObjectCreate(C); 1577 C->data = (void *) (c = PetscNew(Mat_SeqAIJ)); CHKPTRQ(c); 1578 PetscMemcpy(&C->ops,&A->ops,sizeof(struct _MatOps)); 1579 C->destroy = MatDestroy_SeqAIJ; 1580 C->view = MatView_SeqAIJ; 1581 C->factor = A->factor; 1582 c->row = 0; 1583 c->col = 0; 1584 c->indexshift = shift; 1585 C->assembled = PETSC_TRUE; 1586 1587 c->m = C->m = a->m; 1588 c->n = C->n = a->n; 1589 C->M = a->m; 1590 C->N = a->n; 1591 1592 c->imax = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(c->imax); 1593 c->ilen = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(c->ilen); 1594 for ( i=0; i<m; i++ ) { 1595 c->imax[i] = a->imax[i]; 1596 c->ilen[i] = a->ilen[i]; 1597 } 1598 1599 /* allocate the matrix space */ 1600 c->singlemalloc = 1; 1601 len = (m+1)*sizeof(int)+(a->i[m])*(sizeof(Scalar)+sizeof(int)); 1602 c->a = (Scalar *) PetscMalloc( len ); CHKPTRQ(c->a); 1603 c->j = (int *) (c->a + a->i[m] + shift); 1604 c->i = c->j + a->i[m] + shift; 1605 PetscMemcpy(c->i,a->i,(m+1)*sizeof(int)); 1606 if (m > 0) { 1607 PetscMemcpy(c->j,a->j,(a->i[m]+shift)*sizeof(int)); 1608 if (cpvalues == COPY_VALUES) { 1609 PetscMemcpy(c->a,a->a,(a->i[m]+shift)*sizeof(Scalar)); 1610 } 1611 } 1612 1613 PLogObjectMemory(C,len+2*(m+1)*sizeof(int)+sizeof(struct _Mat)+sizeof(Mat_SeqAIJ)); 1614 c->sorted = a->sorted; 1615 c->roworiented = a->roworiented; 1616 c->nonew = a->nonew; 1617 c->ilu_preserve_row_sums = a->ilu_preserve_row_sums; 1618 1619 if (a->diag) { 1620 c->diag = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(c->diag); 1621 PLogObjectMemory(C,(m+1)*sizeof(int)); 1622 for ( i=0; i<m; i++ ) { 1623 c->diag[i] = a->diag[i]; 1624 } 1625 } 1626 else c->diag = 0; 1627 c->inode.limit = a->inode.limit; 1628 c->inode.max_limit = a->inode.max_limit; 1629 if (a->inode.size){ 1630 c->inode.size = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(c->inode.size); 1631 c->inode.node_count = a->inode.node_count; 1632 PetscMemcpy( c->inode.size, a->inode.size, (m+1)*sizeof(int)); 1633 } else { 1634 c->inode.size = 0; 1635 c->inode.node_count = 0; 1636 } 1637 c->nz = a->nz; 1638 c->maxnz = a->maxnz; 1639 c->solve_work = 0; 1640 c->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 1641 1642 *B = C; 1643 return 0; 1644 } 1645 1646 int MatLoad_SeqAIJ(Viewer viewer,MatType type,Mat *A) 1647 { 1648 Mat_SeqAIJ *a; 1649 Mat B; 1650 int i, nz, ierr, fd, header[4],size,*rowlengths = 0,M,N,shift; 1651 MPI_Comm comm; 1652 1653 PetscObjectGetComm((PetscObject) viewer,&comm); 1654 MPI_Comm_size(comm,&size); 1655 if (size > 1) SETERRQ(1,"MatLoad_SeqAIJ:view must have one processor"); 1656 ierr = ViewerBinaryGetDescriptor(viewer,&fd); CHKERRQ(ierr); 1657 ierr = PetscBinaryRead(fd,header,4,BINARY_INT); CHKERRQ(ierr); 1658 if (header[0] != MAT_COOKIE) SETERRQ(1,"MatLoad_SeqAIJ:not matrix object in file"); 1659 M = header[1]; N = header[2]; nz = header[3]; 1660 1661 /* read in row lengths */ 1662 rowlengths = (int*) PetscMalloc( M*sizeof(int) ); CHKPTRQ(rowlengths); 1663 ierr = PetscBinaryRead(fd,rowlengths,M,BINARY_INT); CHKERRQ(ierr); 1664 1665 /* create our matrix */ 1666 ierr = MatCreateSeqAIJ(comm,M,N,0,rowlengths,A); CHKERRQ(ierr); 1667 B = *A; 1668 a = (Mat_SeqAIJ *) B->data; 1669 shift = a->indexshift; 1670 1671 /* read in column indices and adjust for Fortran indexing*/ 1672 ierr = PetscBinaryRead(fd,a->j,nz,BINARY_INT); CHKERRQ(ierr); 1673 if (shift) { 1674 for ( i=0; i<nz; i++ ) { 1675 a->j[i] += 1; 1676 } 1677 } 1678 1679 /* read in nonzero values */ 1680 ierr = PetscBinaryRead(fd,a->a,nz,BINARY_SCALAR); CHKERRQ(ierr); 1681 1682 /* set matrix "i" values */ 1683 a->i[0] = -shift; 1684 for ( i=1; i<= M; i++ ) { 1685 a->i[i] = a->i[i-1] + rowlengths[i-1]; 1686 a->ilen[i-1] = rowlengths[i-1]; 1687 } 1688 PetscFree(rowlengths); 1689 1690 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1691 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1692 return 0; 1693 } 1694 1695 static int MatEqual_SeqAIJ(Mat A,Mat B, PetscTruth* flg) 1696 { 1697 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)B->data; 1698 1699 if (B->type !=MATSEQAIJ)SETERRQ(1,"MatEqual_SeqAIJ:Matrices must be same type"); 1700 1701 /* If the matrix dimensions are not equal, or no of nonzeros or shift */ 1702 if ((a->m != b->m ) || (a->n !=b->n) ||( a->nz != b->nz)|| 1703 (a->indexshift != b->indexshift)) { 1704 *flg = PETSC_FALSE; return 0; 1705 } 1706 1707 /* if the a->i are the same */ 1708 if (PetscMemcmp(a->i,b->i, (a->n+1)*sizeof(int))) { 1709 *flg = PETSC_FALSE; return 0; 1710 } 1711 1712 /* if a->j are the same */ 1713 if (PetscMemcmp(a->j, b->j, (a->nz)*sizeof(int))) { 1714 *flg = PETSC_FALSE; return 0; 1715 } 1716 1717 /* if a->a are the same */ 1718 if (PetscMemcmp(a->a, b->a, (a->nz)*sizeof(Scalar))) { 1719 *flg = PETSC_FALSE; return 0; 1720 } 1721 *flg = PETSC_TRUE; 1722 return 0; 1723 1724 } 1725