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