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