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