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