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