1 2 3 #ifndef lint 4 static char vcid[] = "$Id: aij.c,v 1.190 1996/10/16 21:04:23 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 VecGetArray(ll,&l); VecGetSize(ll,&m); 1092 if (m != a->m) SETERRQ(1,"MatDiagonalScale_SeqAIJ:Left scaling vector wrong length"); 1093 v = a->a; 1094 for ( i=0; i<m; i++ ) { 1095 x = l[i]; 1096 M = a->i[i+1] - a->i[i]; 1097 for ( j=0; j<M; j++ ) { (*v++) *= x;} 1098 } 1099 PLogFlops(nz); 1100 } 1101 if (rr) { 1102 VecGetArray(rr,&r); VecGetSize(rr,&n); 1103 if (n != a->n) SETERRQ(1,"MatDiagonalScale_SeqAIJ:Right scaling vector wrong length"); 1104 v = a->a; jj = a->j; 1105 for ( i=0; i<nz; i++ ) { 1106 (*v++) *= r[*jj++ + shift]; 1107 } 1108 PLogFlops(nz); 1109 } 1110 return 0; 1111 } 1112 1113 static int MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,MatGetSubMatrixCall scall,Mat *B) 1114 { 1115 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data,*c; 1116 int nznew, *smap, i, k, kstart, kend, ierr, oldcols = a->n,*lens; 1117 int row,mat_i,*mat_j,tcol,first,step,*mat_ilen; 1118 register int sum,lensi; 1119 int *irow, *icol, nrows, ncols, shift = a->indexshift,*ssmap; 1120 int *starts,*j_new,*i_new,*aj = a->j, *ai = a->i,ii,*ailen = a->ilen; 1121 Scalar *a_new,*mat_a; 1122 Mat C; 1123 1124 ierr = ISSorted(isrow,(PetscTruth*)&i); 1125 if (!i) SETERRQ(1,"MatGetSubmatrices_SeqAIJ:ISrow is not sorted"); 1126 ierr = ISSorted(iscol,(PetscTruth*)&i); 1127 if (!i) SETERRQ(1,"MatGetSubmatrices_SeqAIJ:IScol is not sorted"); 1128 1129 ierr = ISGetIndices(isrow,&irow); CHKERRQ(ierr); 1130 ierr = ISGetSize(isrow,&nrows); CHKERRQ(ierr); 1131 ierr = ISGetSize(iscol,&ncols); CHKERRQ(ierr); 1132 1133 if (ISStrideGetInfo(iscol,&first,&step) && step == 1) { /* no need to sort */ 1134 /* special case of contiguous rows */ 1135 lens = (int *) PetscMalloc((ncols+nrows+1)*sizeof(int)); CHKPTRQ(lens); 1136 starts = lens + ncols; 1137 /* loop over new rows determining lens and starting points */ 1138 for (i=0; i<nrows; i++) { 1139 kstart = ai[irow[i]]+shift; 1140 kend = kstart + ailen[irow[i]]; 1141 for ( k=kstart; k<kend; k++ ) { 1142 if (aj[k]+shift >= first) { 1143 starts[i] = k; 1144 break; 1145 } 1146 } 1147 sum = 0; 1148 while (k < kend) { 1149 if (aj[k++]+shift >= first+ncols) break; 1150 sum++; 1151 } 1152 lens[i] = sum; 1153 } 1154 /* create submatrix */ 1155 if (scall == MAT_REUSE_MATRIX) { 1156 int n_cols,n_rows; 1157 ierr = MatGetSize(*B,&n_rows,&n_cols); CHKERRQ(ierr); 1158 if (n_rows != nrows || n_cols != ncols) SETERRQ(1,"MatGetSubMatrix_SeqAIJ:"); 1159 ierr = MatZeroEntries(*B); CHKERRQ(ierr); 1160 C = *B; 1161 } 1162 else { 1163 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1164 } 1165 c = (Mat_SeqAIJ*) C->data; 1166 1167 /* loop over rows inserting into submatrix */ 1168 a_new = c->a; 1169 j_new = c->j; 1170 i_new = c->i; 1171 i_new[0] = -shift; 1172 for (i=0; i<nrows; i++) { 1173 ii = starts[i]; 1174 lensi = lens[i]; 1175 for ( k=0; k<lensi; k++ ) { 1176 *j_new++ = aj[ii+k] - first; 1177 } 1178 PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(Scalar)); 1179 a_new += lensi; 1180 i_new[i+1] = i_new[i] + lensi; 1181 c->ilen[i] = lensi; 1182 } 1183 PetscFree(lens); 1184 } 1185 else { 1186 ierr = ISGetIndices(iscol,&icol); CHKERRQ(ierr); 1187 smap = (int *) PetscMalloc((1+oldcols)*sizeof(int)); CHKPTRQ(smap); 1188 ssmap = smap + shift; 1189 lens = (int *) PetscMalloc((1+nrows)*sizeof(int)); CHKPTRQ(lens); 1190 PetscMemzero(smap,oldcols*sizeof(int)); 1191 for ( i=0; i<ncols; i++ ) smap[icol[i]] = i+1; 1192 /* determine lens of each row */ 1193 for (i=0; i<nrows; i++) { 1194 kstart = ai[irow[i]]+shift; 1195 kend = kstart + a->ilen[irow[i]]; 1196 lens[i] = 0; 1197 for ( k=kstart; k<kend; k++ ) { 1198 if (ssmap[aj[k]]) { 1199 lens[i]++; 1200 } 1201 } 1202 } 1203 /* Create and fill new matrix */ 1204 if (scall == MAT_REUSE_MATRIX) { 1205 c = (Mat_SeqAIJ *)((*B)->data); 1206 1207 if (c->m != nrows || c->n != ncols) SETERRQ(1,"MatGetSubMatrix_SeqAIJ:"); 1208 if (PetscMemcmp(c->ilen,lens, c->m *sizeof(int))) { 1209 SETERRQ(1,"MatGetSubmatrices_SeqAIJ:Cannot reuse matrix. wrong no of nonzeros"); 1210 } 1211 PetscMemzero(c->ilen,c->m*sizeof(int)); 1212 C = *B; 1213 } 1214 else { 1215 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1216 } 1217 c = (Mat_SeqAIJ *)(C->data); 1218 for (i=0; i<nrows; i++) { 1219 row = irow[i]; 1220 nznew = 0; 1221 kstart = ai[row]+shift; 1222 kend = kstart + a->ilen[row]; 1223 mat_i = c->i[i]+shift; 1224 mat_j = c->j + mat_i; 1225 mat_a = c->a + mat_i; 1226 mat_ilen = c->ilen + i; 1227 for ( k=kstart; k<kend; k++ ) { 1228 if ((tcol=ssmap[a->j[k]])) { 1229 *mat_j++ = tcol - (!shift); 1230 *mat_a++ = a->a[k]; 1231 (*mat_ilen)++; 1232 1233 } 1234 } 1235 } 1236 /* Free work space */ 1237 ierr = ISRestoreIndices(iscol,&icol); CHKERRQ(ierr); 1238 PetscFree(smap); PetscFree(lens); 1239 } 1240 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1241 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1242 1243 ierr = ISRestoreIndices(isrow,&irow); CHKERRQ(ierr); 1244 *B = C; 1245 return 0; 1246 } 1247 1248 /* 1249 note: This can only work for identity for row and col. It would 1250 be good to check this and otherwise generate an error. 1251 */ 1252 static int MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,double efill,int fill) 1253 { 1254 Mat_SeqAIJ *a = (Mat_SeqAIJ *) inA->data; 1255 int ierr; 1256 Mat outA; 1257 1258 if (fill != 0) SETERRQ(1,"MatILUFactor_SeqAIJ:Only fill=0 supported"); 1259 1260 outA = inA; 1261 inA->factor = FACTOR_LU; 1262 a->row = row; 1263 a->col = col; 1264 1265 a->solve_work = (Scalar *) PetscMalloc( (a->m+1)*sizeof(Scalar)); CHKPTRQ(a->solve_work); 1266 1267 if (!a->diag) { 1268 ierr = MatMarkDiag_SeqAIJ(inA); CHKERRQ(ierr); 1269 } 1270 ierr = MatLUFactorNumeric_SeqAIJ(inA,&outA); CHKERRQ(ierr); 1271 return 0; 1272 } 1273 1274 #include "pinclude/plapack.h" 1275 static int MatScale_SeqAIJ(Scalar *alpha,Mat inA) 1276 { 1277 Mat_SeqAIJ *a = (Mat_SeqAIJ *) inA->data; 1278 int one = 1; 1279 BLscal_( &a->nz, alpha, a->a, &one ); 1280 PLogFlops(a->nz); 1281 return 0; 1282 } 1283 1284 static int MatGetSubMatrices_SeqAIJ(Mat A,int n, IS *irow,IS *icol,MatGetSubMatrixCall scall, 1285 Mat **B) 1286 { 1287 int ierr,i; 1288 1289 if (scall == MAT_INITIAL_MATRIX) { 1290 *B = (Mat *) PetscMalloc( (n+1)*sizeof(Mat) ); CHKPTRQ(*B); 1291 } 1292 1293 for ( i=0; i<n; i++ ) { 1294 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],scall,&(*B)[i]);CHKERRQ(ierr); 1295 } 1296 return 0; 1297 } 1298 1299 static int MatGetBlockSize_SeqAIJ(Mat A, int *bs) 1300 { 1301 *bs = 1; 1302 return 0; 1303 } 1304 1305 static int MatIncreaseOverlap_SeqAIJ(Mat A, int is_max, IS *is, int ov) 1306 { 1307 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1308 int shift, row, i,j,k,l,m,n, *idx,ierr, *nidx, isz, val; 1309 int start, end, *ai, *aj; 1310 char *table; 1311 shift = a->indexshift; 1312 m = a->m; 1313 ai = a->i; 1314 aj = a->j+shift; 1315 1316 if (ov < 0) SETERRQ(1,"MatIncreaseOverlap_SeqAIJ: illegal overlap value used"); 1317 1318 table = (char *) PetscMalloc((m/BITSPERBYTE +1)*sizeof(char)); CHKPTRQ(table); 1319 nidx = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(nidx); 1320 1321 for ( i=0; i<is_max; i++ ) { 1322 /* Initialise the two local arrays */ 1323 isz = 0; 1324 PetscMemzero(table,(m/BITSPERBYTE +1)*sizeof(char)); 1325 1326 /* Extract the indices, assume there can be duplicate entries */ 1327 ierr = ISGetIndices(is[i],&idx); CHKERRQ(ierr); 1328 ierr = ISGetSize(is[i],&n); CHKERRQ(ierr); 1329 1330 /* Enter these into the temp arrays i.e mark table[row], enter row into new index */ 1331 for ( j=0; j<n ; ++j){ 1332 if(!BT_LOOKUP(table, idx[j])) { nidx[isz++] = idx[j];} 1333 } 1334 ierr = ISRestoreIndices(is[i],&idx); CHKERRQ(ierr); 1335 ierr = ISDestroy(is[i]); CHKERRQ(ierr); 1336 1337 k = 0; 1338 for ( j=0; j<ov; j++){ /* for each overlap*/ 1339 n = isz; 1340 for ( ; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 1341 row = nidx[k]; 1342 start = ai[row]; 1343 end = ai[row+1]; 1344 for ( l = start; l<end ; l++){ 1345 val = aj[l] + shift; 1346 if (!BT_LOOKUP(table,val)) {nidx[isz++] = val;} 1347 } 1348 } 1349 } 1350 ierr = ISCreateGeneral(MPI_COMM_SELF, isz, nidx, (is+i)); CHKERRQ(ierr); 1351 } 1352 PetscFree(table); 1353 PetscFree(nidx); 1354 return 0; 1355 } 1356 1357 int MatPrintHelp_SeqAIJ(Mat A) 1358 { 1359 static int called = 0; 1360 MPI_Comm comm = A->comm; 1361 1362 if (called) return 0; else called = 1; 1363 PetscPrintf(comm," Options for MATSEQAIJ and MATMPIAIJ matrix formats (the defaults):\n"); 1364 PetscPrintf(comm," -mat_lu_pivotthreshold <threshold>\n"); 1365 PetscPrintf(comm," -mat_aij_oneindex - internal indices begin at 1 instead of the default 0.\n"); 1366 PetscPrintf(comm," -mat_aij_no_inode - Do not use inodes\n"); 1367 PetscPrintf(comm," -mat_aij_inode_limit <limit> - Set inode limit (max limit=5)\n"); 1368 #if defined(HAVE_ESSL) 1369 PetscPrintf(comm," -mat_aij_essl - Use IBM sparse LU factorization and solve.\n"); 1370 #endif 1371 return 0; 1372 } 1373 static int MatEqual_SeqAIJ(Mat A,Mat B, PetscTruth* flg); 1374 extern int MatFDColoringCreate_SeqAIJ(Mat,ISColoring,MatFDColoring); 1375 extern int MatColoringPatch_SeqAIJ(Mat,int,int *,ISColoring *); 1376 1377 /* -------------------------------------------------------------------*/ 1378 static struct _MatOps MatOps = {MatSetValues_SeqAIJ, 1379 MatGetRow_SeqAIJ,MatRestoreRow_SeqAIJ, 1380 MatMult_SeqAIJ,MatMultAdd_SeqAIJ, 1381 MatMultTrans_SeqAIJ,MatMultTransAdd_SeqAIJ, 1382 MatSolve_SeqAIJ,MatSolveAdd_SeqAIJ, 1383 MatSolveTrans_SeqAIJ,MatSolveTransAdd_SeqAIJ, 1384 MatLUFactor_SeqAIJ,0, 1385 MatRelax_SeqAIJ, 1386 MatTranspose_SeqAIJ, 1387 MatGetInfo_SeqAIJ,MatEqual_SeqAIJ, 1388 MatGetDiagonal_SeqAIJ,MatDiagonalScale_SeqAIJ,MatNorm_SeqAIJ, 1389 0,MatAssemblyEnd_SeqAIJ, 1390 MatCompress_SeqAIJ, 1391 MatSetOption_SeqAIJ,MatZeroEntries_SeqAIJ,MatZeroRows_SeqAIJ, 1392 MatLUFactorSymbolic_SeqAIJ,MatLUFactorNumeric_SeqAIJ,0,0, 1393 MatGetSize_SeqAIJ,MatGetSize_SeqAIJ,MatGetOwnershipRange_SeqAIJ, 1394 MatILUFactorSymbolic_SeqAIJ,0, 1395 0,0,MatConvert_SeqAIJ, 1396 MatConvertSameType_SeqAIJ,0,0, 1397 MatILUFactor_SeqAIJ,0,0, 1398 MatGetSubMatrices_SeqAIJ,MatIncreaseOverlap_SeqAIJ, 1399 MatGetValues_SeqAIJ,0, 1400 MatPrintHelp_SeqAIJ, 1401 MatScale_SeqAIJ,0,0, 1402 MatILUDTFactor_SeqAIJ, 1403 MatGetBlockSize_SeqAIJ, 1404 MatGetRowIJ_SeqAIJ, 1405 MatRestoreRowIJ_SeqAIJ, 1406 MatGetColumnIJ_SeqAIJ, 1407 MatRestoreColumnIJ_SeqAIJ, 1408 MatFDColoringCreate_SeqAIJ, 1409 MatColoringPatch_SeqAIJ}; 1410 1411 extern int MatUseSuperLU_SeqAIJ(Mat); 1412 extern int MatUseEssl_SeqAIJ(Mat); 1413 extern int MatUseDXML_SeqAIJ(Mat); 1414 1415 /*@C 1416 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 1417 (the default parallel PETSc format). For good matrix assembly performance 1418 the user should preallocate the matrix storage by setting the parameter nz 1419 (or the array nzz). By setting these parameters accurately, performance 1420 during matrix assembly can be increased by more than a factor of 50. 1421 1422 Input Parameters: 1423 . comm - MPI communicator, set to MPI_COMM_SELF 1424 . m - number of rows 1425 . n - number of columns 1426 . nz - number of nonzeros per row (same for all rows) 1427 . nzz - array containing the number of nonzeros in the various rows 1428 (possibly different for each row) or PETSC_NULL 1429 1430 Output Parameter: 1431 . A - the matrix 1432 1433 Notes: 1434 The AIJ format (also called the Yale sparse matrix format or 1435 compressed row storage), is fully compatible with standard Fortran 77 1436 storage. That is, the stored row and column indices can begin at 1437 either one (as in Fortran) or zero. See the users' manual for details. 1438 1439 Specify the preallocated storage with either nz or nnz (not both). 1440 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1441 allocation. For large problems you MUST preallocate memory or you 1442 will get TERRIBLE performance, see the users' manual chapter on 1443 matrices and the file $(PETSC_DIR)/Performance. 1444 1445 By default, this format uses inodes (identical nodes) when possible, to 1446 improve numerical efficiency of Matrix vector products and solves. We 1447 search for consecutive rows with the same nonzero structure, thereby 1448 reusing matrix information to achieve increased efficiency. 1449 1450 Options Database Keys: 1451 $ -mat_aij_no_inode - Do not use inodes 1452 $ -mat_aij_inode_limit <limit> - Set inode limit. 1453 $ (max limit=5) 1454 $ -mat_aij_oneindex - Internally use indexing starting at 1 1455 $ rather than 0. Note: When calling MatSetValues(), 1456 $ the user still MUST index entries starting at 0! 1457 1458 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues() 1459 @*/ 1460 int MatCreateSeqAIJ(MPI_Comm comm,int m,int n,int nz,int *nnz, Mat *A) 1461 { 1462 Mat B; 1463 Mat_SeqAIJ *b; 1464 int i, len, ierr, flg,size; 1465 1466 MPI_Comm_size(comm,&size); 1467 if (size > 1) SETERRQ(1,"MatCreateSeqAIJ:Comm must be of size 1"); 1468 1469 *A = 0; 1470 PetscHeaderCreate(B,_Mat,MAT_COOKIE,MATSEQAIJ,comm); 1471 PLogObjectCreate(B); 1472 B->data = (void *) (b = PetscNew(Mat_SeqAIJ)); CHKPTRQ(b); 1473 PetscMemzero(b,sizeof(Mat_SeqAIJ)); 1474 PetscMemcpy(&B->ops,&MatOps,sizeof(struct _MatOps)); 1475 B->destroy = MatDestroy_SeqAIJ; 1476 B->view = MatView_SeqAIJ; 1477 B->factor = 0; 1478 B->lupivotthreshold = 1.0; 1479 ierr = OptionsGetDouble(PETSC_NULL,"-mat_lu_pivotthreshold",&B->lupivotthreshold, 1480 &flg); CHKERRQ(ierr); 1481 b->ilu_preserve_row_sums = PETSC_FALSE; 1482 ierr = OptionsHasName(PETSC_NULL,"-pc_ilu_preserve_row_sums", 1483 (int*) &b->ilu_preserve_row_sums); CHKERRQ(ierr); 1484 b->row = 0; 1485 b->col = 0; 1486 b->indexshift = 0; 1487 b->reallocs = 0; 1488 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_oneindex", &flg); CHKERRQ(ierr); 1489 if (flg) b->indexshift = -1; 1490 1491 b->m = m; B->m = m; B->M = m; 1492 b->n = n; B->n = n; B->N = n; 1493 b->imax = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(b->imax); 1494 if (nnz == PETSC_NULL) { 1495 if (nz == PETSC_DEFAULT) nz = 10; 1496 else if (nz <= 0) nz = 1; 1497 for ( i=0; i<m; i++ ) b->imax[i] = nz; 1498 nz = nz*m; 1499 } 1500 else { 1501 nz = 0; 1502 for ( i=0; i<m; i++ ) {b->imax[i] = nnz[i]; nz += nnz[i];} 1503 } 1504 1505 /* allocate the matrix space */ 1506 len = nz*(sizeof(int) + sizeof(Scalar)) + (b->m+1)*sizeof(int); 1507 b->a = (Scalar *) PetscMalloc( len ); CHKPTRQ(b->a); 1508 b->j = (int *) (b->a + nz); 1509 PetscMemzero(b->j,nz*sizeof(int)); 1510 b->i = b->j + nz; 1511 b->singlemalloc = 1; 1512 1513 b->i[0] = -b->indexshift; 1514 for (i=1; i<m+1; i++) { 1515 b->i[i] = b->i[i-1] + b->imax[i-1]; 1516 } 1517 1518 /* b->ilen will count nonzeros in each row so far. */ 1519 b->ilen = (int *) PetscMalloc((m+1)*sizeof(int)); 1520 PLogObjectMemory(B,len+2*(m+1)*sizeof(int)+sizeof(struct _Mat)+sizeof(Mat_SeqAIJ)); 1521 for ( i=0; i<b->m; i++ ) { b->ilen[i] = 0;} 1522 1523 b->nz = 0; 1524 b->maxnz = nz; 1525 b->sorted = 0; 1526 b->roworiented = 1; 1527 b->nonew = 0; 1528 b->diag = 0; 1529 b->solve_work = 0; 1530 b->spptr = 0; 1531 b->inode.node_count = 0; 1532 b->inode.size = 0; 1533 b->inode.limit = 5; 1534 b->inode.max_limit = 5; 1535 B->info.nz_unneeded = (double)b->maxnz; 1536 1537 *A = B; 1538 1539 /* SuperLU is not currently supported through PETSc */ 1540 #if defined(HAVE_SUPERLU) 1541 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_superlu", &flg); CHKERRQ(ierr); 1542 if (flg) { ierr = MatUseSuperLU_SeqAIJ(B); CHKERRQ(ierr); } 1543 #endif 1544 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_essl", &flg); CHKERRQ(ierr); 1545 if (flg) { ierr = MatUseEssl_SeqAIJ(B); CHKERRQ(ierr); } 1546 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_dxml", &flg); CHKERRQ(ierr); 1547 if (flg) { 1548 if (!b->indexshift) SETERRQ(1,"MatCreateSeqAIJ:need -mat_aij_oneindex with -mat_aij_dxml"); 1549 ierr = MatUseDXML_SeqAIJ(B); CHKERRQ(ierr); 1550 } 1551 ierr = OptionsHasName(PETSC_NULL,"-help", &flg); CHKERRQ(ierr); 1552 if (flg) {ierr = MatPrintHelp(B); CHKERRQ(ierr); } 1553 return 0; 1554 } 1555 1556 int MatConvertSameType_SeqAIJ(Mat A,Mat *B,int cpvalues) 1557 { 1558 Mat C; 1559 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ *) A->data; 1560 int i,len, m = a->m,shift = a->indexshift; 1561 1562 *B = 0; 1563 PetscHeaderCreate(C,_Mat,MAT_COOKIE,MATSEQAIJ,A->comm); 1564 PLogObjectCreate(C); 1565 C->data = (void *) (c = PetscNew(Mat_SeqAIJ)); CHKPTRQ(c); 1566 PetscMemcpy(&C->ops,&A->ops,sizeof(struct _MatOps)); 1567 C->destroy = MatDestroy_SeqAIJ; 1568 C->view = MatView_SeqAIJ; 1569 C->factor = A->factor; 1570 c->row = 0; 1571 c->col = 0; 1572 c->indexshift = shift; 1573 C->assembled = PETSC_TRUE; 1574 1575 c->m = C->m = a->m; 1576 c->n = C->n = a->n; 1577 C->M = a->m; 1578 C->N = a->n; 1579 1580 c->imax = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(c->imax); 1581 c->ilen = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(c->ilen); 1582 for ( i=0; i<m; i++ ) { 1583 c->imax[i] = a->imax[i]; 1584 c->ilen[i] = a->ilen[i]; 1585 } 1586 1587 /* allocate the matrix space */ 1588 c->singlemalloc = 1; 1589 len = (m+1)*sizeof(int)+(a->i[m])*(sizeof(Scalar)+sizeof(int)); 1590 c->a = (Scalar *) PetscMalloc( len ); CHKPTRQ(c->a); 1591 c->j = (int *) (c->a + a->i[m] + shift); 1592 c->i = c->j + a->i[m] + shift; 1593 PetscMemcpy(c->i,a->i,(m+1)*sizeof(int)); 1594 if (m > 0) { 1595 PetscMemcpy(c->j,a->j,(a->i[m]+shift)*sizeof(int)); 1596 if (cpvalues == COPY_VALUES) { 1597 PetscMemcpy(c->a,a->a,(a->i[m]+shift)*sizeof(Scalar)); 1598 } 1599 } 1600 1601 PLogObjectMemory(C,len+2*(m+1)*sizeof(int)+sizeof(struct _Mat)+sizeof(Mat_SeqAIJ)); 1602 c->sorted = a->sorted; 1603 c->roworiented = a->roworiented; 1604 c->nonew = a->nonew; 1605 c->ilu_preserve_row_sums = a->ilu_preserve_row_sums; 1606 1607 if (a->diag) { 1608 c->diag = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(c->diag); 1609 PLogObjectMemory(C,(m+1)*sizeof(int)); 1610 for ( i=0; i<m; i++ ) { 1611 c->diag[i] = a->diag[i]; 1612 } 1613 } 1614 else c->diag = 0; 1615 c->inode.limit = a->inode.limit; 1616 c->inode.max_limit = a->inode.max_limit; 1617 if (a->inode.size){ 1618 c->inode.size = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(c->inode.size); 1619 c->inode.node_count = a->inode.node_count; 1620 PetscMemcpy( c->inode.size, a->inode.size, (m+1)*sizeof(int)); 1621 } else { 1622 c->inode.size = 0; 1623 c->inode.node_count = 0; 1624 } 1625 c->nz = a->nz; 1626 c->maxnz = a->maxnz; 1627 c->solve_work = 0; 1628 c->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 1629 1630 *B = C; 1631 return 0; 1632 } 1633 1634 int MatLoad_SeqAIJ(Viewer viewer,MatType type,Mat *A) 1635 { 1636 Mat_SeqAIJ *a; 1637 Mat B; 1638 int i, nz, ierr, fd, header[4],size,*rowlengths = 0,M,N,shift; 1639 MPI_Comm comm; 1640 1641 PetscObjectGetComm((PetscObject) viewer,&comm); 1642 MPI_Comm_size(comm,&size); 1643 if (size > 1) SETERRQ(1,"MatLoad_SeqAIJ:view must have one processor"); 1644 ierr = ViewerBinaryGetDescriptor(viewer,&fd); CHKERRQ(ierr); 1645 ierr = PetscBinaryRead(fd,header,4,BINARY_INT); CHKERRQ(ierr); 1646 if (header[0] != MAT_COOKIE) SETERRQ(1,"MatLoad_SeqAIJ:not matrix object in file"); 1647 M = header[1]; N = header[2]; nz = header[3]; 1648 1649 /* read in row lengths */ 1650 rowlengths = (int*) PetscMalloc( M*sizeof(int) ); CHKPTRQ(rowlengths); 1651 ierr = PetscBinaryRead(fd,rowlengths,M,BINARY_INT); CHKERRQ(ierr); 1652 1653 /* create our matrix */ 1654 ierr = MatCreateSeqAIJ(comm,M,N,0,rowlengths,A); CHKERRQ(ierr); 1655 B = *A; 1656 a = (Mat_SeqAIJ *) B->data; 1657 shift = a->indexshift; 1658 1659 /* read in column indices and adjust for Fortran indexing*/ 1660 ierr = PetscBinaryRead(fd,a->j,nz,BINARY_INT); CHKERRQ(ierr); 1661 if (shift) { 1662 for ( i=0; i<nz; i++ ) { 1663 a->j[i] += 1; 1664 } 1665 } 1666 1667 /* read in nonzero values */ 1668 ierr = PetscBinaryRead(fd,a->a,nz,BINARY_SCALAR); CHKERRQ(ierr); 1669 1670 /* set matrix "i" values */ 1671 a->i[0] = -shift; 1672 for ( i=1; i<= M; i++ ) { 1673 a->i[i] = a->i[i-1] + rowlengths[i-1]; 1674 a->ilen[i-1] = rowlengths[i-1]; 1675 } 1676 PetscFree(rowlengths); 1677 1678 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1679 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1680 return 0; 1681 } 1682 1683 static int MatEqual_SeqAIJ(Mat A,Mat B, PetscTruth* flg) 1684 { 1685 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)B->data; 1686 1687 if (B->type !=MATSEQAIJ)SETERRQ(1,"MatEqual_SeqAIJ:Matrices must be same type"); 1688 1689 /* If the matrix dimensions are not equal, or no of nonzeros or shift */ 1690 if ((a->m != b->m ) || (a->n !=b->n) ||( a->nz != b->nz)|| 1691 (a->indexshift != b->indexshift)) { 1692 *flg = PETSC_FALSE; return 0; 1693 } 1694 1695 /* if the a->i are the same */ 1696 if (PetscMemcmp(a->i,b->i, (a->n+1)*sizeof(int))) { 1697 *flg = PETSC_FALSE; return 0; 1698 } 1699 1700 /* if a->j are the same */ 1701 if (PetscMemcmp(a->j, b->j, (a->nz)*sizeof(int))) { 1702 *flg = PETSC_FALSE; return 0; 1703 } 1704 1705 /* if a->a are the same */ 1706 if (PetscMemcmp(a->a, b->a, (a->nz)*sizeof(Scalar))) { 1707 *flg = PETSC_FALSE; return 0; 1708 } 1709 *flg = PETSC_TRUE; 1710 return 0; 1711 1712 } 1713