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