1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: aij.c,v 1.273 1998/07/13 18:52:55 balay Exp bsmith $"; 3 #endif 4 5 /* 6 Defines the basic matrix operations for the AIJ (compressed row) 7 matrix storage format. 8 */ 9 10 #include "pinclude/pviewer.h" 11 #include "sys.h" 12 #include "src/mat/impls/aij/seq/aij.h" 13 #include "src/vec/vecimpl.h" 14 #include "src/inline/spops.h" 15 #include "src/inline/dot.h" 16 #include "src/inline/bitarray.h" 17 18 /* 19 Basic AIJ format ILU based on drop tolerance 20 */ 21 #undef __FUNC__ 22 #define __FUNC__ "MatILUDTFactor_SeqAIJ" 23 int MatILUDTFactor_SeqAIJ(Mat A,double dt,int maxnz,IS row,IS col,Mat *fact) 24 { 25 /* Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; */ 26 int ierr = 1; 27 28 PetscFunctionBegin; 29 SETERRQ(ierr,0,"Not implemented"); 30 #if !defined(USE_PETSC_DEBUG) 31 PetscFunctionReturn(0); 32 #endif 33 } 34 35 extern int MatToSymmetricIJ_SeqAIJ(int,int*,int*,int,int,int**,int**); 36 37 #undef __FUNC__ 38 #define __FUNC__ "MatGetRowIJ_SeqAIJ" 39 int MatGetRowIJ_SeqAIJ(Mat A,int oshift,PetscTruth symmetric,int *m,int **ia,int **ja, 40 PetscTruth *done) 41 { 42 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 43 int ierr,i,ishift; 44 45 PetscFunctionBegin; 46 *m = A->m; 47 if (!ia) PetscFunctionReturn(0); 48 ishift = a->indexshift; 49 if (symmetric) { 50 ierr = MatToSymmetricIJ_SeqAIJ(a->m,a->i,a->j,ishift,oshift,ia,ja); CHKERRQ(ierr); 51 } else if (oshift == 0 && ishift == -1) { 52 int nz = a->i[a->m]; 53 /* malloc space and subtract 1 from i and j indices */ 54 *ia = (int *) PetscMalloc( (a->m+1)*sizeof(int) ); CHKPTRQ(*ia); 55 *ja = (int *) PetscMalloc( (nz+1)*sizeof(int) ); CHKPTRQ(*ja); 56 for ( i=0; i<nz; i++ ) (*ja)[i] = a->j[i] - 1; 57 for ( i=0; i<a->m+1; i++ ) (*ia)[i] = a->i[i] - 1; 58 } else if (oshift == 1 && ishift == 0) { 59 int nz = a->i[a->m] + 1; 60 /* malloc space and add 1 to i and j indices */ 61 *ia = (int *) PetscMalloc( (a->m+1)*sizeof(int) ); CHKPTRQ(*ia); 62 *ja = (int *) PetscMalloc( (nz+1)*sizeof(int) ); CHKPTRQ(*ja); 63 for ( i=0; i<nz; i++ ) (*ja)[i] = a->j[i] + 1; 64 for ( i=0; i<a->m+1; i++ ) (*ia)[i] = a->i[i] + 1; 65 } else { 66 *ia = a->i; *ja = a->j; 67 } 68 69 PetscFunctionReturn(0); 70 } 71 72 #undef __FUNC__ 73 #define __FUNC__ "MatRestoreRowIJ_SeqAIJ" 74 int MatRestoreRowIJ_SeqAIJ(Mat A,int oshift,PetscTruth symmetric,int *n,int **ia,int **ja, 75 PetscTruth *done) 76 { 77 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 78 int ishift = a->indexshift; 79 80 PetscFunctionBegin; 81 if (!ia) PetscFunctionReturn(0); 82 if (symmetric || (oshift == 0 && ishift == -1) || (oshift == 1 && ishift == 0)) { 83 PetscFree(*ia); 84 PetscFree(*ja); 85 } 86 PetscFunctionReturn(0); 87 } 88 89 #undef __FUNC__ 90 #define __FUNC__ "MatGetColumnIJ_SeqAIJ" 91 int MatGetColumnIJ_SeqAIJ(Mat A,int oshift,PetscTruth symmetric,int *nn,int **ia,int **ja, 92 PetscTruth *done) 93 { 94 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 95 int ierr,i,ishift = a->indexshift,*collengths,*cia,*cja,n = A->n,m = A->m; 96 int nz = a->i[m]+ishift,row,*jj,mr,col; 97 98 PetscFunctionBegin; 99 *nn = A->n; 100 if (!ia) PetscFunctionReturn(0); 101 if (symmetric) { 102 ierr = MatToSymmetricIJ_SeqAIJ(a->m,a->i,a->j,ishift,oshift,ia,ja); CHKERRQ(ierr); 103 } else { 104 collengths = (int *) PetscMalloc( (n+1)*sizeof(int) ); CHKPTRQ(collengths); 105 PetscMemzero(collengths,n*sizeof(int)); 106 cia = (int *) PetscMalloc( (n+1)*sizeof(int) ); CHKPTRQ(cia); 107 cja = (int *) PetscMalloc( (nz+1)*sizeof(int) ); CHKPTRQ(cja); 108 jj = a->j; 109 for ( i=0; i<nz; i++ ) { 110 collengths[jj[i] + ishift]++; 111 } 112 cia[0] = oshift; 113 for ( i=0; i<n; i++) { 114 cia[i+1] = cia[i] + collengths[i]; 115 } 116 PetscMemzero(collengths,n*sizeof(int)); 117 jj = a->j; 118 for ( row=0; row<m; row++ ) { 119 mr = a->i[row+1] - a->i[row]; 120 for ( i=0; i<mr; i++ ) { 121 col = *jj++ + ishift; 122 cja[cia[col] + collengths[col]++ - oshift] = row + oshift; 123 } 124 } 125 PetscFree(collengths); 126 *ia = cia; *ja = cja; 127 } 128 129 PetscFunctionReturn(0); 130 } 131 132 #undef __FUNC__ 133 #define __FUNC__ "MatRestoreColumnIJ_SeqAIJ" 134 int MatRestoreColumnIJ_SeqAIJ(Mat A,int oshift,PetscTruth symmetric,int *n,int **ia, 135 int **ja,PetscTruth *done) 136 { 137 PetscFunctionBegin; 138 if (!ia) PetscFunctionReturn(0); 139 140 PetscFree(*ia); 141 PetscFree(*ja); 142 143 PetscFunctionReturn(0); 144 } 145 146 #define CHUNKSIZE 15 147 148 #undef __FUNC__ 149 #define __FUNC__ "MatSetValues_SeqAIJ" 150 int MatSetValues_SeqAIJ(Mat A,int m,int *im,int n,int *in,Scalar *v,InsertMode is) 151 { 152 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 153 int *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax, N, sorted = a->sorted; 154 int *imax = a->imax, *ai = a->i, *ailen = a->ilen,roworiented = a->roworiented; 155 int *aj = a->j, nonew = a->nonew,shift = a->indexshift; 156 Scalar *ap,value, *aa = a->a; 157 158 PetscFunctionBegin; 159 for ( k=0; k<m; k++ ) { /* loop over added rows */ 160 row = im[k]; 161 #if defined(USE_PETSC_BOPT_g) 162 if (row < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative row"); 163 if (row >= a->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row too large"); 164 #endif 165 rp = aj + ai[row] + shift; ap = aa + ai[row] + shift; 166 rmax = imax[row]; nrow = ailen[row]; 167 low = 0; 168 for ( l=0; l<n; l++ ) { /* loop over added columns */ 169 #if defined(USE_PETSC_BOPT_g) 170 if (in[l] < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative column"); 171 if (in[l] >= a->n) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Column too large"); 172 #endif 173 col = in[l] - shift; 174 if (roworiented) { 175 value = *v++; 176 } else { 177 value = v[k + l*m]; 178 } 179 if (!sorted) low = 0; high = nrow; 180 while (high-low > 5) { 181 t = (low+high)/2; 182 if (rp[t] > col) high = t; 183 else low = t; 184 } 185 for ( i=low; i<high; i++ ) { 186 if (rp[i] > col) break; 187 if (rp[i] == col) { 188 if (is == ADD_VALUES) ap[i] += value; 189 else ap[i] = value; 190 goto noinsert; 191 } 192 } 193 if (nonew == 1) goto noinsert; 194 else if (nonew == -1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Inserting a new nonzero in the matrix"); 195 if (nrow >= rmax) { 196 /* there is no extra room in row, therefore enlarge */ 197 int new_nz = ai[a->m] + CHUNKSIZE,len,*new_i,*new_j; 198 Scalar *new_a; 199 200 if (nonew == -2) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Inserting a new nonzero in the matrix"); 201 202 /* malloc new storage space */ 203 len = new_nz*(sizeof(int)+sizeof(Scalar))+(a->m+1)*sizeof(int); 204 new_a = (Scalar *) PetscMalloc( len ); CHKPTRQ(new_a); 205 new_j = (int *) (new_a + new_nz); 206 new_i = new_j + new_nz; 207 208 /* copy over old data into new slots */ 209 for ( ii=0; ii<row+1; ii++ ) {new_i[ii] = ai[ii];} 210 for ( ii=row+1; ii<a->m+1; ii++ ) {new_i[ii] = ai[ii]+CHUNKSIZE;} 211 PetscMemcpy(new_j,aj,(ai[row]+nrow+shift)*sizeof(int)); 212 len = (new_nz - CHUNKSIZE - ai[row] - nrow - shift); 213 PetscMemcpy(new_j+ai[row]+shift+nrow+CHUNKSIZE,aj+ai[row]+shift+nrow, 214 len*sizeof(int)); 215 PetscMemcpy(new_a,aa,(ai[row]+nrow+shift)*sizeof(Scalar)); 216 PetscMemcpy(new_a+ai[row]+shift+nrow+CHUNKSIZE,aa+ai[row]+shift+nrow, 217 len*sizeof(Scalar)); 218 /* free up old matrix storage */ 219 PetscFree(a->a); 220 if (!a->singlemalloc) {PetscFree(a->i);PetscFree(a->j);} 221 aa = a->a = new_a; ai = a->i = new_i; aj = a->j = new_j; 222 a->singlemalloc = 1; 223 224 rp = aj + ai[row] + shift; ap = aa + ai[row] + shift; 225 rmax = imax[row] = imax[row] + CHUNKSIZE; 226 PLogObjectMemory(A,CHUNKSIZE*(sizeof(int) + sizeof(Scalar))); 227 a->maxnz += CHUNKSIZE; 228 a->reallocs++; 229 } 230 N = nrow++ - 1; a->nz++; 231 /* shift up all the later entries in this row */ 232 for ( ii=N; ii>=i; ii-- ) { 233 rp[ii+1] = rp[ii]; 234 ap[ii+1] = ap[ii]; 235 } 236 rp[i] = col; 237 ap[i] = value; 238 noinsert:; 239 low = i + 1; 240 } 241 ailen[row] = nrow; 242 } 243 PetscFunctionReturn(0); 244 } 245 246 #undef __FUNC__ 247 #define __FUNC__ "MatGetValues_SeqAIJ" 248 int MatGetValues_SeqAIJ(Mat A,int m,int *im,int n,int *in,Scalar *v) 249 { 250 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 251 int *rp, k, low, high, t, row, nrow, i, col, l, *aj = a->j; 252 int *ai = a->i, *ailen = a->ilen, shift = a->indexshift; 253 Scalar *ap, *aa = a->a, zero = 0.0; 254 255 PetscFunctionBegin; 256 for ( k=0; k<m; k++ ) { /* loop over rows */ 257 row = im[k]; 258 if (row < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative row"); 259 if (row >= a->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row too large"); 260 rp = aj + ai[row] + shift; ap = aa + ai[row] + shift; 261 nrow = ailen[row]; 262 for ( l=0; l<n; l++ ) { /* loop over columns */ 263 if (in[l] < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Negative column"); 264 if (in[l] >= a->n) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Column too large"); 265 col = in[l] - shift; 266 high = nrow; low = 0; /* assume unsorted */ 267 while (high-low > 5) { 268 t = (low+high)/2; 269 if (rp[t] > col) high = t; 270 else low = t; 271 } 272 for ( i=low; i<high; i++ ) { 273 if (rp[i] > col) break; 274 if (rp[i] == col) { 275 *v++ = ap[i]; 276 goto finished; 277 } 278 } 279 *v++ = zero; 280 finished:; 281 } 282 } 283 PetscFunctionReturn(0); 284 } 285 286 287 #undef __FUNC__ 288 #define __FUNC__ "MatView_SeqAIJ_Binary" 289 int MatView_SeqAIJ_Binary(Mat A,Viewer viewer) 290 { 291 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 292 int i, fd, *col_lens, ierr; 293 294 PetscFunctionBegin; 295 ierr = ViewerBinaryGetDescriptor(viewer,&fd); CHKERRQ(ierr); 296 col_lens = (int *) PetscMalloc( (4+a->m)*sizeof(int) ); CHKPTRQ(col_lens); 297 col_lens[0] = MAT_COOKIE; 298 col_lens[1] = a->m; 299 col_lens[2] = a->n; 300 col_lens[3] = a->nz; 301 302 /* store lengths of each row and write (including header) to file */ 303 for ( i=0; i<a->m; i++ ) { 304 col_lens[4+i] = a->i[i+1] - a->i[i]; 305 } 306 ierr = PetscBinaryWrite(fd,col_lens,4+a->m,PETSC_INT,1); CHKERRQ(ierr); 307 PetscFree(col_lens); 308 309 /* store column indices (zero start index) */ 310 if (a->indexshift) { 311 for ( i=0; i<a->nz; i++ ) a->j[i]--; 312 } 313 ierr = PetscBinaryWrite(fd,a->j,a->nz,PETSC_INT,0); CHKERRQ(ierr); 314 if (a->indexshift) { 315 for ( i=0; i<a->nz; i++ ) a->j[i]++; 316 } 317 318 /* store nonzero values */ 319 ierr = PetscBinaryWrite(fd,a->a,a->nz,PETSC_SCALAR,0); CHKERRQ(ierr); 320 PetscFunctionReturn(0); 321 } 322 323 #undef __FUNC__ 324 #define __FUNC__ "MatView_SeqAIJ_ASCII" 325 int MatView_SeqAIJ_ASCII(Mat A,Viewer viewer) 326 { 327 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 328 int ierr, i,j, m = a->m, shift = a->indexshift, format, flg1,flg2; 329 FILE *fd; 330 char *outputname; 331 332 PetscFunctionBegin; 333 ierr = ViewerASCIIGetPointer(viewer,&fd); CHKERRQ(ierr); 334 ierr = ViewerFileGetOutputname_Private(viewer,&outputname); CHKERRQ(ierr); 335 ierr = ViewerGetFormat(viewer,&format); 336 if (format == VIEWER_FORMAT_ASCII_INFO) { 337 PetscFunctionReturn(0); 338 } else if (format == VIEWER_FORMAT_ASCII_INFO_LONG) { 339 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_no_inode",&flg1); CHKERRQ(ierr); 340 ierr = OptionsHasName(PETSC_NULL,"-mat_no_unroll",&flg2); CHKERRQ(ierr); 341 if (flg1 || flg2) fprintf(fd," not using I-node routines\n"); 342 else fprintf(fd," using I-node routines: found %d nodes, limit used is %d\n", 343 a->inode.node_count,a->inode.limit); 344 } else if (format == VIEWER_FORMAT_ASCII_MATLAB) { 345 int nofinalvalue = 0; 346 if ((a->i[m] == a->i[m-1]) || (a->j[a->nz-1] != a->n-!shift)) { 347 nofinalvalue = 1; 348 } 349 fprintf(fd,"%% Size = %d %d \n",m,a->n); 350 fprintf(fd,"%% Nonzeros = %d \n",a->nz); 351 fprintf(fd,"zzz = zeros(%d,3);\n",a->nz+nofinalvalue); 352 fprintf(fd,"zzz = [\n"); 353 354 for (i=0; i<m; i++) { 355 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 356 #if defined(USE_PETSC_COMPLEX) 357 fprintf(fd,"%d %d %18.16e + %18.16e i \n",i+1,a->j[j]+!shift,PetscReal(a->a[j]),PetscImaginary(a->a[j])); 358 #else 359 fprintf(fd,"%d %d %18.16e\n", i+1, a->j[j]+!shift, a->a[j]); 360 #endif 361 } 362 } 363 if (nofinalvalue) { 364 fprintf(fd,"%d %d %18.16e\n", m, a->n, 0.0); 365 } 366 fprintf(fd,"];\n %s = spconvert(zzz);\n",outputname); 367 } else if (format == VIEWER_FORMAT_ASCII_COMMON) { 368 for ( i=0; i<m; i++ ) { 369 fprintf(fd,"row %d:",i); 370 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 371 #if defined(USE_PETSC_COMPLEX) 372 if (PetscImaginary(a->a[j]) > 0.0 && PetscReal(a->a[j]) != 0.0) 373 fprintf(fd," %d %g + %g i",a->j[j]+shift,PetscReal(a->a[j]),PetscImaginary(a->a[j])); 374 else if (PetscImaginary(a->a[j]) < 0.0 && PetscReal(a->a[j]) != 0.0) 375 fprintf(fd," %d %g - %g i",a->j[j]+shift,PetscReal(a->a[j]),-PetscImaginary(a->a[j])); 376 else if (PetscReal(a->a[j]) != 0.0) 377 fprintf(fd," %d %g ",a->j[j]+shift,PetscReal(a->a[j])); 378 #else 379 if (a->a[j] != 0.0) fprintf(fd," %d %g ",a->j[j]+shift,a->a[j]); 380 #endif 381 } 382 fprintf(fd,"\n"); 383 } 384 } 385 else if (format == VIEWER_FORMAT_ASCII_SYMMODU) { 386 int nzd=0, fshift=1, *sptr; 387 sptr = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(sptr); 388 for ( i=0; i<m; i++ ) { 389 sptr[i] = nzd+1; 390 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 391 if (a->j[j] >= i) { 392 #if defined(USE_PETSC_COMPLEX) 393 if (PetscImaginary(a->a[j]) != 0.0 || PetscReal(a->a[j]) != 0.0) nzd++; 394 #else 395 if (a->a[j] != 0.0) nzd++; 396 #endif 397 } 398 } 399 } 400 sptr[m] = nzd+1; 401 fprintf(fd," %d %d\n\n",m,nzd); 402 for ( i=0; i<m+1; i+=6 ) { 403 if (i+4<m) fprintf(fd," %d %d %d %d %d %d\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3],sptr[i+4],sptr[i+5]); 404 else if (i+3<m) fprintf(fd," %d %d %d %d %d\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3],sptr[i+4]); 405 else if (i+2<m) fprintf(fd," %d %d %d %d\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3]); 406 else if (i+1<m) fprintf(fd," %d %d %d\n",sptr[i],sptr[i+1],sptr[i+2]); 407 else if (i<m) fprintf(fd," %d %d\n",sptr[i],sptr[i+1]); 408 else fprintf(fd," %d\n",sptr[i]); 409 } 410 fprintf(fd,"\n"); 411 PetscFree(sptr); 412 for ( i=0; i<m; i++ ) { 413 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 414 if (a->j[j] >= i) fprintf(fd," %d ",a->j[j]+fshift); 415 } 416 fprintf(fd,"\n"); 417 } 418 fprintf(fd,"\n"); 419 for ( i=0; i<m; i++ ) { 420 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 421 if (a->j[j] >= i) { 422 #if defined(USE_PETSC_COMPLEX) 423 if (PetscImaginary(a->a[j]) != 0.0 || PetscReal(a->a[j]) != 0.0) 424 fprintf(fd," %18.16e %18.16e ",PetscReal(a->a[j]),PetscImaginary(a->a[j])); 425 #else 426 if (a->a[j] != 0.0) fprintf(fd," %18.16e ",a->a[j]); 427 #endif 428 } 429 } 430 fprintf(fd,"\n"); 431 } 432 } else { 433 for ( i=0; i<m; i++ ) { 434 fprintf(fd,"row %d:",i); 435 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 436 #if defined(USE_PETSC_COMPLEX) 437 if (PetscImaginary(a->a[j]) > 0.0) { 438 fprintf(fd," %d %g + %g i",a->j[j]+shift,PetscReal(a->a[j]),PetscImaginary(a->a[j])); 439 } else if (PetscImaginary(a->a[j]) < 0.0) { 440 fprintf(fd," %d %g - %g i",a->j[j]+shift,PetscReal(a->a[j]),-PetscImaginary(a->a[j])); 441 } else { 442 fprintf(fd," %d %g ",a->j[j]+shift,PetscReal(a->a[j])); 443 } 444 #else 445 fprintf(fd," %d %g ",a->j[j]+shift,a->a[j]); 446 #endif 447 } 448 fprintf(fd,"\n"); 449 } 450 } 451 fflush(fd); 452 PetscFunctionReturn(0); 453 } 454 455 #undef __FUNC__ 456 #define __FUNC__ "MatView_SeqAIJ_Draw_Zoom" 457 int MatView_SeqAIJ_Draw_Zoom(Draw draw,void *Aa) 458 { 459 Mat A = (Mat) Aa; 460 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 461 int ierr, i,j, m = a->m, shift = a->indexshift,color; 462 int format; 463 double xl,yl,xr,yr,x_l,x_r,y_l,y_r,maxv = 0.0; 464 Viewer viewer; 465 466 PetscFunctionBegin; 467 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*) &viewer);CHKERRQ(ierr); 468 ierr = ViewerGetFormat(viewer,&format); CHKERRQ(ierr); 469 470 ierr = DrawGetCoordinates(draw,&xl,&yl,&xr,&yr); CHKERRQ(ierr); 471 /* loop over matrix elements drawing boxes */ 472 473 if (format != VIEWER_FORMAT_DRAW_CONTOUR) { 474 /* Blue for negative, Cyan for zero and Red for positive */ 475 color = DRAW_BLUE; 476 for ( i=0; i<m; i++ ) { 477 y_l = m - i - 1.0; y_r = y_l + 1.0; 478 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 479 x_l = a->j[j] + shift; x_r = x_l + 1.0; 480 #if defined(USE_PETSC_COMPLEX) 481 if (PetscReal(a->a[j]) >= 0.) continue; 482 #else 483 if (a->a[j] >= 0.) continue; 484 #endif 485 ierr = DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 486 } 487 } 488 color = DRAW_CYAN; 489 for ( i=0; i<m; i++ ) { 490 y_l = m - i - 1.0; y_r = y_l + 1.0; 491 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 492 x_l = a->j[j] + shift; x_r = x_l + 1.0; 493 if (a->a[j] != 0.) continue; 494 ierr = DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 495 } 496 } 497 color = DRAW_RED; 498 for ( i=0; i<m; i++ ) { 499 y_l = m - i - 1.0; y_r = y_l + 1.0; 500 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 501 x_l = a->j[j] + shift; x_r = x_l + 1.0; 502 #if defined(USE_PETSC_COMPLEX) 503 if (PetscReal(a->a[j]) <= 0.) continue; 504 #else 505 if (a->a[j] <= 0.) continue; 506 #endif 507 ierr = DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 508 } 509 } 510 } else { 511 /* use contour shading to indicate magnitude of values */ 512 /* first determine max of all nonzero values */ 513 int nz = a->nz,count; 514 Draw popup; 515 double scale; 516 517 for ( i=0; i<nz; i++ ) { 518 if (PetscAbsScalar(a->a[i]) > maxv) maxv = PetscAbsScalar(a->a[i]); 519 } 520 scale = (245.0 - DRAW_BASIC_COLORS)/maxv; 521 ierr = DrawGetPopup(draw,&popup); CHKERRQ(ierr); 522 ierr = DrawScalePopup(popup,0.0,maxv); CHKERRQ(ierr); 523 count = 0; 524 for ( i=0; i<m; i++ ) { 525 y_l = m - i - 1.0; y_r = y_l + 1.0; 526 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 527 x_l = a->j[j] + shift; x_r = x_l + 1.0; 528 color = DRAW_BASIC_COLORS + (int)(scale*PetscAbsScalar(a->a[count])); 529 ierr = DrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 530 count++; 531 } 532 } 533 } 534 PetscFunctionReturn(0); 535 } 536 537 #undef __FUNC__ 538 #define __FUNC__ "MatView_SeqAIJ_Draw" 539 int MatView_SeqAIJ_Draw(Mat A,Viewer viewer) 540 { 541 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 542 int ierr; 543 Draw draw; 544 double xr,yr,xl,yl,h,w; 545 546 PetscTruth isnull; 547 548 PetscFunctionBegin; 549 ierr = ViewerDrawGetDraw(viewer,&draw); CHKERRQ(ierr); 550 ierr = DrawIsNull(draw,&isnull); CHKERRQ(ierr); 551 if (isnull) PetscFunctionReturn(0); 552 553 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 554 xr = a->n; yr = a->m; h = yr/10.0; w = xr/10.0; 555 xr += w; yr += h; xl = -w; yl = -h; 556 ierr = DrawSetCoordinates(draw,xl,yl,xr,yr); CHKERRQ(ierr); 557 ierr = DrawZoom(draw,MatView_SeqAIJ_Draw_Zoom,A); CHKERRQ(ierr); 558 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);CHKERRQ(ierr); 559 PetscFunctionReturn(0); 560 } 561 562 #undef __FUNC__ 563 #define __FUNC__ "MatView_SeqAIJ" 564 int MatView_SeqAIJ(Mat A,Viewer viewer) 565 { 566 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 567 ViewerType vtype; 568 int ierr; 569 570 PetscFunctionBegin; 571 ierr = ViewerGetType(viewer,&vtype); CHKERRQ(ierr); 572 if (vtype == MATLAB_VIEWER) { 573 ierr = ViewerMatlabPutSparse_Private(viewer,a->m,a->n,a->nz,a->a,a->i,a->j); CHKERRQ(ierr); 574 } else if (vtype == ASCII_FILE_VIEWER || vtype == ASCII_FILES_VIEWER){ 575 ierr = MatView_SeqAIJ_ASCII(A,viewer); CHKERRQ(ierr); 576 } else if (vtype == BINARY_FILE_VIEWER) { 577 ierr = MatView_SeqAIJ_Binary(A,viewer); CHKERRQ(ierr); 578 } else if (vtype == DRAW_VIEWER) { 579 ierr = MatView_SeqAIJ_Draw(A,viewer); CHKERRQ(ierr); 580 } else { 581 SETERRQ(1,1,"Viewer type not supported by PETSc object"); 582 } 583 PetscFunctionReturn(0); 584 } 585 586 extern int Mat_AIJ_CheckInode(Mat); 587 #undef __FUNC__ 588 #define __FUNC__ "MatAssemblyEnd_SeqAIJ" 589 int MatAssemblyEnd_SeqAIJ(Mat A,MatAssemblyType mode) 590 { 591 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 592 int fshift = 0,i,j,*ai = a->i, *aj = a->j, *imax = a->imax,ierr; 593 int m = a->m, *ip, N, *ailen = a->ilen,shift = a->indexshift,rmax = 0; 594 Scalar *aa = a->a, *ap; 595 596 PetscFunctionBegin; 597 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 598 599 if (m) rmax = ailen[0]; /* determine row with most nonzeros */ 600 for ( i=1; i<m; i++ ) { 601 /* move each row back by the amount of empty slots (fshift) before it*/ 602 fshift += imax[i-1] - ailen[i-1]; 603 rmax = PetscMax(rmax,ailen[i]); 604 if (fshift) { 605 ip = aj + ai[i] + shift; ap = aa + ai[i] + shift; 606 N = ailen[i]; 607 for ( j=0; j<N; j++ ) { 608 ip[j-fshift] = ip[j]; 609 ap[j-fshift] = ap[j]; 610 } 611 } 612 ai[i] = ai[i-1] + ailen[i-1]; 613 } 614 if (m) { 615 fshift += imax[m-1] - ailen[m-1]; 616 ai[m] = ai[m-1] + ailen[m-1]; 617 } 618 /* reset ilen and imax for each row */ 619 for ( i=0; i<m; i++ ) { 620 ailen[i] = imax[i] = ai[i+1] - ai[i]; 621 } 622 a->nz = ai[m] + shift; 623 624 /* diagonals may have moved, so kill the diagonal pointers */ 625 if (fshift && a->diag) { 626 PetscFree(a->diag); 627 PLogObjectMemory(A,-(m+1)*sizeof(int)); 628 a->diag = 0; 629 } 630 PLogInfo(A,"MatAssemblyEnd_SeqAIJ:Matrix size: %d X %d; storage space: %d unneeded, %d used\n", 631 m,a->n,fshift,a->nz); 632 PLogInfo(A,"MatAssemblyEnd_SeqAIJ:Number of mallocs during MatSetValues is %d\n", 633 a->reallocs); 634 PLogInfo(A,"MatAssemblyEnd_SeqAIJ:Most nonzeros in any row is %d\n",rmax); 635 a->reallocs = 0; 636 A->info.nz_unneeded = (double)fshift; 637 638 /* check out for identical nodes. If found, use inode functions */ 639 ierr = Mat_AIJ_CheckInode(A); CHKERRQ(ierr); 640 PetscFunctionReturn(0); 641 } 642 643 #undef __FUNC__ 644 #define __FUNC__ "MatZeroEntries_SeqAIJ" 645 int MatZeroEntries_SeqAIJ(Mat A) 646 { 647 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 648 649 PetscFunctionBegin; 650 PetscMemzero(a->a,(a->i[a->m]+a->indexshift)*sizeof(Scalar)); 651 PetscFunctionReturn(0); 652 } 653 654 #undef __FUNC__ 655 #define __FUNC__ "MatDestroy_SeqAIJ" 656 int MatDestroy_SeqAIJ(Mat A) 657 { 658 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 659 int ierr; 660 661 PetscFunctionBegin; 662 if (--A->refct > 0) PetscFunctionReturn(0); 663 664 if (A->mapping) { 665 ierr = ISLocalToGlobalMappingDestroy(A->mapping); CHKERRQ(ierr); 666 } 667 if (A->bmapping) { 668 ierr = ISLocalToGlobalMappingDestroy(A->bmapping); CHKERRQ(ierr); 669 } 670 #if defined(USE_PETSC_LOG) 671 PLogObjectState((PetscObject)A,"Rows=%d, Cols=%d, NZ=%d",a->m,a->n,a->nz); 672 #endif 673 PetscFree(a->a); 674 if (!a->singlemalloc) { PetscFree(a->i); PetscFree(a->j);} 675 if (a->diag) PetscFree(a->diag); 676 if (a->ilen) PetscFree(a->ilen); 677 if (a->imax) PetscFree(a->imax); 678 if (a->solve_work) PetscFree(a->solve_work); 679 if (a->inode.size) PetscFree(a->inode.size); 680 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} 681 PetscFree(a); 682 683 PLogObjectDestroy(A); 684 PetscHeaderDestroy(A); 685 PetscFunctionReturn(0); 686 } 687 688 #undef __FUNC__ 689 #define __FUNC__ "MatCompress_SeqAIJ" 690 int MatCompress_SeqAIJ(Mat A) 691 { 692 PetscFunctionBegin; 693 PetscFunctionReturn(0); 694 } 695 696 #undef __FUNC__ 697 #define __FUNC__ "MatSetOption_SeqAIJ" 698 int MatSetOption_SeqAIJ(Mat A,MatOption op) 699 { 700 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 701 702 PetscFunctionBegin; 703 if (op == MAT_ROW_ORIENTED) a->roworiented = 1; 704 else if (op == MAT_COLUMN_ORIENTED) a->roworiented = 0; 705 else if (op == MAT_COLUMNS_SORTED) a->sorted = 1; 706 else if (op == MAT_COLUMNS_UNSORTED) a->sorted = 0; 707 else if (op == MAT_NO_NEW_NONZERO_LOCATIONS) a->nonew = 1; 708 else if (op == MAT_NEW_NONZERO_LOCATION_ERROR) a->nonew = -1; 709 else if (op == MAT_NEW_NONZERO_ALLOCATION_ERROR) a->nonew = -2; 710 else if (op == MAT_YES_NEW_NONZERO_LOCATIONS) a->nonew = 0; 711 else if (op == MAT_ROWS_SORTED || 712 op == MAT_ROWS_UNSORTED || 713 op == MAT_SYMMETRIC || 714 op == MAT_STRUCTURALLY_SYMMETRIC || 715 op == MAT_YES_NEW_DIAGONALS || 716 op == MAT_IGNORE_OFF_PROC_ENTRIES|| 717 op == MAT_USE_HASH_TABLE) 718 PLogInfo(A,"MatSetOption_SeqAIJ:Option ignored\n"); 719 else if (op == MAT_NO_NEW_DIAGONALS) { 720 SETERRQ(PETSC_ERR_SUP,0,"MAT_NO_NEW_DIAGONALS"); 721 } else if (op == MAT_INODE_LIMIT_1) a->inode.limit = 1; 722 else if (op == MAT_INODE_LIMIT_2) a->inode.limit = 2; 723 else if (op == MAT_INODE_LIMIT_3) a->inode.limit = 3; 724 else if (op == MAT_INODE_LIMIT_4) a->inode.limit = 4; 725 else if (op == MAT_INODE_LIMIT_5) a->inode.limit = 5; 726 else SETERRQ(PETSC_ERR_SUP,0,"unknown option"); 727 PetscFunctionReturn(0); 728 } 729 730 #undef __FUNC__ 731 #define __FUNC__ "MatGetDiagonal_SeqAIJ" 732 int MatGetDiagonal_SeqAIJ(Mat A,Vec v) 733 { 734 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 735 int i,j, n,shift = a->indexshift,ierr; 736 Scalar *x, zero = 0.0; 737 738 PetscFunctionBegin; 739 ierr = VecSet(&zero,v);CHKERRQ(ierr); 740 ierr = VecGetArray(v,&x);;CHKERRQ(ierr); 741 ierr = VecGetLocalSize(v,&n);;CHKERRQ(ierr); 742 if (n != a->m) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Nonconforming matrix and vector"); 743 for ( i=0; i<a->m; i++ ) { 744 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 745 if (a->j[j]+shift == i) { 746 x[i] = a->a[j]; 747 break; 748 } 749 } 750 } 751 ierr = VecRestoreArray(v,&x);;CHKERRQ(ierr); 752 PetscFunctionReturn(0); 753 } 754 755 /* -------------------------------------------------------*/ 756 /* Should check that shapes of vectors and matrices match */ 757 /* -------------------------------------------------------*/ 758 #undef __FUNC__ 759 #define __FUNC__ "MatMultTrans_SeqAIJ" 760 int MatMultTrans_SeqAIJ(Mat A,Vec xx,Vec yy) 761 { 762 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 763 Scalar *x, *y, *v, alpha; 764 int ierr,m = a->m, n, i, *idx, shift = a->indexshift; 765 766 PetscFunctionBegin; 767 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 768 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 769 PetscMemzero(y,a->n*sizeof(Scalar)); 770 y = y + shift; /* shift for Fortran start by 1 indexing */ 771 for ( i=0; i<m; i++ ) { 772 idx = a->j + a->i[i] + shift; 773 v = a->a + a->i[i] + shift; 774 n = a->i[i+1] - a->i[i]; 775 alpha = x[i]; 776 while (n-->0) {y[*idx++] += alpha * *v++;} 777 } 778 PLogFlops(2*a->nz - a->n); 779 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 780 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 781 PetscFunctionReturn(0); 782 } 783 784 #undef __FUNC__ 785 #define __FUNC__ "MatMultTransAdd_SeqAIJ" 786 int MatMultTransAdd_SeqAIJ(Mat A,Vec xx,Vec zz,Vec yy) 787 { 788 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 789 Scalar *x, *y, *v, alpha; 790 int ierr,m = a->m, n, i, *idx,shift = a->indexshift; 791 792 PetscFunctionBegin; 793 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 794 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 795 if (zz != yy) VecCopy(zz,yy); 796 y = y + shift; /* shift for Fortran start by 1 indexing */ 797 for ( i=0; i<m; i++ ) { 798 idx = a->j + a->i[i] + shift; 799 v = a->a + a->i[i] + shift; 800 n = a->i[i+1] - a->i[i]; 801 alpha = x[i]; 802 while (n-->0) {y[*idx++] += alpha * *v++;} 803 } 804 PLogFlops(2*a->nz); 805 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 806 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 807 PetscFunctionReturn(0); 808 } 809 810 #undef __FUNC__ 811 #define __FUNC__ "MatMult_SeqAIJ" 812 int MatMult_SeqAIJ(Mat A,Vec xx,Vec yy) 813 { 814 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 815 Scalar *x, *y, *v, sum; 816 int ierr,m = a->m, *idx, shift = a->indexshift,*ii; 817 #if !defined(USE_FORTRAN_KERNEL_MULTAIJ) 818 int n, i, jrow,j; 819 #endif 820 821 #if defined(HAVE_PRAGMA_DISJOINT) 822 #pragma disjoint(*x,*y,*v) 823 #endif 824 825 PetscFunctionBegin; 826 ierr = VecGetArray(xx,&x); CHKERRQ(ierr); 827 ierr = VecGetArray(yy,&y); CHKERRQ(ierr); 828 x = x + shift; /* shift for Fortran start by 1 indexing */ 829 idx = a->j; 830 v = a->a; 831 ii = a->i; 832 #if defined(USE_FORTRAN_KERNEL_MULTAIJ) 833 fortranmultaij_(&m,x,ii,idx+shift,v+shift,y); 834 #else 835 v += shift; /* shift for Fortran start by 1 indexing */ 836 idx += shift; 837 for ( i=0; i<m; i++ ) { 838 jrow = ii[i]; 839 n = ii[i+1] - jrow; 840 sum = 0.0; 841 for ( j=0; j<n; j++) { 842 sum += v[jrow]*x[idx[jrow]]; jrow++; 843 } 844 y[i] = sum; 845 } 846 #endif 847 PLogFlops(2*a->nz - m); 848 ierr = VecRestoreArray(xx,&x); CHKERRQ(ierr); 849 ierr = VecRestoreArray(yy,&y); CHKERRQ(ierr); 850 PetscFunctionReturn(0); 851 } 852 853 #undef __FUNC__ 854 #define __FUNC__ "MatMultAdd_SeqAIJ" 855 int MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 856 { 857 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 858 Scalar *x, *y, *z, *v, sum; 859 int ierr,m = a->m, *idx, shift = a->indexshift,*ii; 860 #if !defined(USE_FORTRAN_KERNEL_MULTADDAIJ) 861 int n,i,jrow,j; 862 #endif 863 864 PetscFunctionBegin; 865 ierr = VecGetArray(xx,&x); CHKERRQ(ierr); 866 ierr = VecGetArray(yy,&y); CHKERRQ(ierr); 867 ierr = VecGetArray(zz,&z); CHKERRQ(ierr); 868 x = x + shift; /* shift for Fortran start by 1 indexing */ 869 idx = a->j; 870 v = a->a; 871 ii = a->i; 872 #if defined(USE_FORTRAN_KERNEL_MULTADDAIJ) 873 fortranmultaddaij_(&m,x,ii,idx+shift,v+shift,y,z); 874 #else 875 v += shift; /* shift for Fortran start by 1 indexing */ 876 idx += shift; 877 for ( i=0; i<m; i++ ) { 878 jrow = ii[i]; 879 n = ii[i+1] - jrow; 880 sum = y[i]; 881 for ( j=0; j<n; j++) { 882 sum += v[jrow]*x[idx[jrow]]; jrow++; 883 } 884 z[i] = sum; 885 } 886 #endif 887 PLogFlops(2*a->nz); 888 ierr = VecRestoreArray(xx,&x); CHKERRQ(ierr); 889 ierr = VecRestoreArray(yy,&y); CHKERRQ(ierr); 890 ierr = VecRestoreArray(zz,&z); CHKERRQ(ierr); 891 PetscFunctionReturn(0); 892 } 893 894 /* 895 Adds diagonal pointers to sparse matrix structure. 896 */ 897 898 #undef __FUNC__ 899 #define __FUNC__ "MatMarkDiag_SeqAIJ" 900 int MatMarkDiag_SeqAIJ(Mat A) 901 { 902 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 903 int i,j, *diag, m = a->m,shift = a->indexshift; 904 905 PetscFunctionBegin; 906 diag = (int *) PetscMalloc( (m+1)*sizeof(int)); CHKPTRQ(diag); 907 PLogObjectMemory(A,(m+1)*sizeof(int)); 908 for ( i=0; i<a->m; i++ ) { 909 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 910 if (a->j[j]+shift == i) { 911 diag[i] = j - shift; 912 break; 913 } 914 } 915 } 916 a->diag = diag; 917 PetscFunctionReturn(0); 918 } 919 920 #undef __FUNC__ 921 #define __FUNC__ "MatRelax_SeqAIJ" 922 int MatRelax_SeqAIJ(Mat A,Vec bb,double omega,MatSORType flag, 923 double fshift,int its,Vec xx) 924 { 925 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 926 Scalar *x, *b, *bs, d, *xs, sum, *v = a->a,*t,scale,*ts, *xb; 927 int ierr, *idx, *diag,n = a->n, m = a->m, i, shift = a->indexshift; 928 929 PetscFunctionBegin; 930 ierr = VecGetArray(xx,&x); CHKERRQ(ierr); 931 ierr = VecGetArray(bb,&b);CHKERRQ(ierr); 932 if (!a->diag) {ierr = MatMarkDiag_SeqAIJ(A);CHKERRQ(ierr);} 933 diag = a->diag; 934 xs = x + shift; /* shifted by one for index start of a or a->j*/ 935 if (flag == SOR_APPLY_UPPER) { 936 /* apply ( U + D/omega) to the vector */ 937 bs = b + shift; 938 for ( i=0; i<m; i++ ) { 939 d = fshift + a->a[diag[i] + shift]; 940 n = a->i[i+1] - diag[i] - 1; 941 idx = a->j + diag[i] + (!shift); 942 v = a->a + diag[i] + (!shift); 943 sum = b[i]*d/omega; 944 SPARSEDENSEDOT(sum,bs,v,idx,n); 945 x[i] = sum; 946 } 947 PetscFunctionReturn(0); 948 } 949 if (flag == SOR_APPLY_LOWER) { 950 SETERRQ(PETSC_ERR_SUP,0,"SOR_APPLY_LOWER is not done"); 951 } else if (flag & SOR_EISENSTAT) { 952 /* Let A = L + U + D; where L is lower trianglar, 953 U is upper triangular, E is diagonal; This routine applies 954 955 (L + E)^{-1} A (U + E)^{-1} 956 957 to a vector efficiently using Eisenstat's trick. This is for 958 the case of SSOR preconditioner, so E is D/omega where omega 959 is the relaxation factor. 960 */ 961 t = (Scalar *) PetscMalloc( m*sizeof(Scalar) ); CHKPTRQ(t); 962 scale = (2.0/omega) - 1.0; 963 964 /* x = (E + U)^{-1} b */ 965 for ( i=m-1; i>=0; i-- ) { 966 d = fshift + a->a[diag[i] + shift]; 967 n = a->i[i+1] - diag[i] - 1; 968 idx = a->j + diag[i] + (!shift); 969 v = a->a + diag[i] + (!shift); 970 sum = b[i]; 971 SPARSEDENSEMDOT(sum,xs,v,idx,n); 972 x[i] = omega*(sum/d); 973 } 974 975 /* t = b - (2*E - D)x */ 976 v = a->a; 977 for ( i=0; i<m; i++ ) { t[i] = b[i] - scale*(v[*diag++ + shift])*x[i]; } 978 979 /* t = (E + L)^{-1}t */ 980 ts = t + shift; /* shifted by one for index start of a or a->j*/ 981 diag = a->diag; 982 for ( i=0; i<m; i++ ) { 983 d = fshift + a->a[diag[i]+shift]; 984 n = diag[i] - a->i[i]; 985 idx = a->j + a->i[i] + shift; 986 v = a->a + a->i[i] + shift; 987 sum = t[i]; 988 SPARSEDENSEMDOT(sum,ts,v,idx,n); 989 t[i] = omega*(sum/d); 990 } 991 992 /* x = x + t */ 993 for ( i=0; i<m; i++ ) { x[i] += t[i]; } 994 PetscFree(t); 995 PetscFunctionReturn(0); 996 } 997 if (flag & SOR_ZERO_INITIAL_GUESS) { 998 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 999 for ( i=0; i<m; i++ ) { 1000 d = fshift + a->a[diag[i]+shift]; 1001 n = diag[i] - a->i[i]; 1002 idx = a->j + a->i[i] + shift; 1003 v = a->a + a->i[i] + shift; 1004 sum = b[i]; 1005 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1006 x[i] = omega*(sum/d); 1007 } 1008 xb = x; 1009 } else xb = b; 1010 if ((flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) && 1011 (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP)) { 1012 for ( i=0; i<m; i++ ) { 1013 x[i] *= a->a[diag[i]+shift]; 1014 } 1015 } 1016 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1017 for ( i=m-1; i>=0; i-- ) { 1018 d = fshift + a->a[diag[i] + shift]; 1019 n = a->i[i+1] - diag[i] - 1; 1020 idx = a->j + diag[i] + (!shift); 1021 v = a->a + diag[i] + (!shift); 1022 sum = xb[i]; 1023 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1024 x[i] = omega*(sum/d); 1025 } 1026 } 1027 its--; 1028 } 1029 while (its--) { 1030 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1031 for ( i=0; i<m; i++ ) { 1032 d = fshift + a->a[diag[i]+shift]; 1033 n = a->i[i+1] - a->i[i]; 1034 idx = a->j + a->i[i] + shift; 1035 v = a->a + a->i[i] + shift; 1036 sum = b[i]; 1037 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1038 x[i] = (1. - omega)*x[i] + omega*(sum + a->a[diag[i]+shift]*x[i])/d; 1039 } 1040 } 1041 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1042 for ( i=m-1; i>=0; i-- ) { 1043 d = fshift + a->a[diag[i] + shift]; 1044 n = a->i[i+1] - a->i[i]; 1045 idx = a->j + a->i[i] + shift; 1046 v = a->a + a->i[i] + shift; 1047 sum = b[i]; 1048 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1049 x[i] = (1. - omega)*x[i] + omega*(sum + a->a[diag[i]+shift]*x[i])/d; 1050 } 1051 } 1052 } 1053 PetscFunctionReturn(0); 1054 } 1055 1056 #undef __FUNC__ 1057 #define __FUNC__ "MatGetInfo_SeqAIJ" 1058 int MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1059 { 1060 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1061 1062 PetscFunctionBegin; 1063 info->rows_global = (double)a->m; 1064 info->columns_global = (double)a->n; 1065 info->rows_local = (double)a->m; 1066 info->columns_local = (double)a->n; 1067 info->block_size = 1.0; 1068 info->nz_allocated = (double)a->maxnz; 1069 info->nz_used = (double)a->nz; 1070 info->nz_unneeded = (double)(a->maxnz - a->nz); 1071 /* if (info->nz_unneeded != A->info.nz_unneeded) 1072 printf("space descrepancy: maxnz-nz = %d, nz_unneeded = %d\n",(int)info->nz_unneeded,(int)A->info.nz_unneeded); */ 1073 info->assemblies = (double)A->num_ass; 1074 info->mallocs = (double)a->reallocs; 1075 info->memory = A->mem; 1076 if (A->factor) { 1077 info->fill_ratio_given = A->info.fill_ratio_given; 1078 info->fill_ratio_needed = A->info.fill_ratio_needed; 1079 info->factor_mallocs = A->info.factor_mallocs; 1080 } else { 1081 info->fill_ratio_given = 0; 1082 info->fill_ratio_needed = 0; 1083 info->factor_mallocs = 0; 1084 } 1085 PetscFunctionReturn(0); 1086 } 1087 1088 extern int MatLUFactorSymbolic_SeqAIJ(Mat,IS,IS,double,Mat*); 1089 extern int MatLUFactorNumeric_SeqAIJ(Mat,Mat*); 1090 extern int MatLUFactor_SeqAIJ(Mat,IS,IS,double); 1091 extern int MatSolve_SeqAIJ(Mat,Vec,Vec); 1092 extern int MatSolveAdd_SeqAIJ(Mat,Vec,Vec,Vec); 1093 extern int MatSolveTrans_SeqAIJ(Mat,Vec,Vec); 1094 extern int MatSolveTransAdd_SeqAIJ(Mat,Vec,Vec,Vec); 1095 1096 #undef __FUNC__ 1097 #define __FUNC__ "MatZeroRows_SeqAIJ" 1098 int MatZeroRows_SeqAIJ(Mat A,IS is,Scalar *diag) 1099 { 1100 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1101 int i,ierr,N, *rows,m = a->m - 1,shift = a->indexshift; 1102 1103 PetscFunctionBegin; 1104 ierr = ISGetSize(is,&N); CHKERRQ(ierr); 1105 ierr = ISGetIndices(is,&rows); CHKERRQ(ierr); 1106 if (diag) { 1107 for ( i=0; i<N; i++ ) { 1108 if (rows[i] < 0 || rows[i] > m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"row out of range"); 1109 if (a->ilen[rows[i]] > 0) { /* in case row was completely empty */ 1110 a->ilen[rows[i]] = 1; 1111 a->a[a->i[rows[i]]+shift] = *diag; 1112 a->j[a->i[rows[i]]+shift] = rows[i]+shift; 1113 } else { 1114 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],diag,INSERT_VALUES);CHKERRQ(ierr); 1115 } 1116 } 1117 } else { 1118 for ( i=0; i<N; i++ ) { 1119 if (rows[i] < 0 || rows[i] > m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"row out of range"); 1120 a->ilen[rows[i]] = 0; 1121 } 1122 } 1123 ISRestoreIndices(is,&rows); 1124 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1125 PetscFunctionReturn(0); 1126 } 1127 1128 #undef __FUNC__ 1129 #define __FUNC__ "MatGetSize_SeqAIJ" 1130 int MatGetSize_SeqAIJ(Mat A,int *m,int *n) 1131 { 1132 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1133 1134 PetscFunctionBegin; 1135 if (m) *m = a->m; 1136 if (n) *n = a->n; 1137 PetscFunctionReturn(0); 1138 } 1139 1140 #undef __FUNC__ 1141 #define __FUNC__ "MatGetOwnershipRange_SeqAIJ" 1142 int MatGetOwnershipRange_SeqAIJ(Mat A,int *m,int *n) 1143 { 1144 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1145 1146 PetscFunctionBegin; 1147 *m = 0; *n = a->m; 1148 PetscFunctionReturn(0); 1149 } 1150 1151 #undef __FUNC__ 1152 #define __FUNC__ "MatGetRow_SeqAIJ" 1153 int MatGetRow_SeqAIJ(Mat A,int row,int *nz,int **idx,Scalar **v) 1154 { 1155 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1156 int *itmp,i,shift = a->indexshift; 1157 1158 PetscFunctionBegin; 1159 if (row < 0 || row >= a->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row out of range"); 1160 1161 *nz = a->i[row+1] - a->i[row]; 1162 if (v) *v = a->a + a->i[row] + shift; 1163 if (idx) { 1164 itmp = a->j + a->i[row] + shift; 1165 if (*nz && shift) { 1166 *idx = (int *) PetscMalloc( (*nz)*sizeof(int) ); CHKPTRQ(*idx); 1167 for ( i=0; i<(*nz); i++ ) {(*idx)[i] = itmp[i] + shift;} 1168 } else if (*nz) { 1169 *idx = itmp; 1170 } 1171 else *idx = 0; 1172 } 1173 PetscFunctionReturn(0); 1174 } 1175 1176 #undef __FUNC__ 1177 #define __FUNC__ "MatRestoreRow_SeqAIJ" 1178 int MatRestoreRow_SeqAIJ(Mat A,int row,int *nz,int **idx,Scalar **v) 1179 { 1180 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1181 1182 PetscFunctionBegin; 1183 if (idx) {if (*idx && a->indexshift) PetscFree(*idx);} 1184 PetscFunctionReturn(0); 1185 } 1186 1187 #undef __FUNC__ 1188 #define __FUNC__ "MatNorm_SeqAIJ" 1189 int MatNorm_SeqAIJ(Mat A,NormType type,double *norm) 1190 { 1191 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1192 Scalar *v = a->a; 1193 double sum = 0.0; 1194 int i, j,shift = a->indexshift; 1195 1196 PetscFunctionBegin; 1197 if (type == NORM_FROBENIUS) { 1198 for (i=0; i<a->nz; i++ ) { 1199 #if defined(USE_PETSC_COMPLEX) 1200 sum += PetscReal(PetscConj(*v)*(*v)); v++; 1201 #else 1202 sum += (*v)*(*v); v++; 1203 #endif 1204 } 1205 *norm = sqrt(sum); 1206 } else if (type == NORM_1) { 1207 double *tmp; 1208 int *jj = a->j; 1209 tmp = (double *) PetscMalloc( (a->n+1)*sizeof(double) ); CHKPTRQ(tmp); 1210 PetscMemzero(tmp,a->n*sizeof(double)); 1211 *norm = 0.0; 1212 for ( j=0; j<a->nz; j++ ) { 1213 tmp[*jj++ + shift] += PetscAbsScalar(*v); v++; 1214 } 1215 for ( j=0; j<a->n; j++ ) { 1216 if (tmp[j] > *norm) *norm = tmp[j]; 1217 } 1218 PetscFree(tmp); 1219 } else if (type == NORM_INFINITY) { 1220 *norm = 0.0; 1221 for ( j=0; j<a->m; j++ ) { 1222 v = a->a + a->i[j] + shift; 1223 sum = 0.0; 1224 for ( i=0; i<a->i[j+1]-a->i[j]; i++ ) { 1225 sum += PetscAbsScalar(*v); v++; 1226 } 1227 if (sum > *norm) *norm = sum; 1228 } 1229 } else { 1230 SETERRQ(PETSC_ERR_SUP,0,"No support for two norm"); 1231 } 1232 PetscFunctionReturn(0); 1233 } 1234 1235 #undef __FUNC__ 1236 #define __FUNC__ "MatTranspose_SeqAIJ" 1237 int MatTranspose_SeqAIJ(Mat A,Mat *B) 1238 { 1239 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1240 Mat C; 1241 int i, ierr, *aj = a->j, *ai = a->i, m = a->m, len, *col; 1242 int shift = a->indexshift; 1243 Scalar *array = a->a; 1244 1245 PetscFunctionBegin; 1246 if (B == PETSC_NULL && m != a->n) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Square matrix only for in-place"); 1247 col = (int *) PetscMalloc((1+a->n)*sizeof(int)); CHKPTRQ(col); 1248 PetscMemzero(col,(1+a->n)*sizeof(int)); 1249 if (shift) { 1250 for ( i=0; i<ai[m]-1; i++ ) aj[i] -= 1; 1251 } 1252 for ( i=0; i<ai[m]+shift; i++ ) col[aj[i]] += 1; 1253 ierr = MatCreateSeqAIJ(A->comm,a->n,m,0,col,&C); CHKERRQ(ierr); 1254 PetscFree(col); 1255 for ( i=0; i<m; i++ ) { 1256 len = ai[i+1]-ai[i]; 1257 ierr = MatSetValues(C,len,aj,1,&i,array,INSERT_VALUES); CHKERRQ(ierr); 1258 array += len; aj += len; 1259 } 1260 if (shift) { 1261 for ( i=0; i<ai[m]-1; i++ ) aj[i] += 1; 1262 } 1263 1264 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1265 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1266 1267 if (B != PETSC_NULL) { 1268 *B = C; 1269 } else { 1270 PetscOps *Abops; 1271 MatOps Aops; 1272 1273 /* This isn't really an in-place transpose */ 1274 PetscFree(a->a); 1275 if (!a->singlemalloc) {PetscFree(a->i); PetscFree(a->j);} 1276 if (a->diag) PetscFree(a->diag); 1277 if (a->ilen) PetscFree(a->ilen); 1278 if (a->imax) PetscFree(a->imax); 1279 if (a->solve_work) PetscFree(a->solve_work); 1280 if (a->inode.size) PetscFree(a->inode.size); 1281 PetscFree(a); 1282 1283 /* 1284 This is horrible, horrible code. We need to keep the 1285 the bops and ops Structures, copy everything from C 1286 including the function pointers.. 1287 */ 1288 PetscMemcpy(A->ops,C->ops,sizeof(struct _MatOps)); 1289 PetscMemcpy(A->bops,C->bops,sizeof(PetscOps)); 1290 Abops = A->bops; 1291 Aops = A->ops; 1292 PetscMemcpy(A,C,sizeof(struct _p_Mat)); 1293 A->bops = Abops; 1294 A->ops = Aops; 1295 A->qlist = 0; 1296 1297 PetscHeaderDestroy(C); 1298 } 1299 PetscFunctionReturn(0); 1300 } 1301 1302 #undef __FUNC__ 1303 #define __FUNC__ "MatDiagonalScale_SeqAIJ" 1304 int MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 1305 { 1306 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1307 Scalar *l,*r,x,*v; 1308 int ierr,i,j,m = a->m, n = a->n, M, nz = a->nz, *jj,shift = a->indexshift; 1309 1310 PetscFunctionBegin; 1311 if (ll) { 1312 /* The local size is used so that VecMPI can be passed to this routine 1313 by MatDiagonalScale_MPIAIJ */ 1314 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 1315 if (m != a->m) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Left scaling vector wrong length"); 1316 ierr = VecGetArray(ll,&l); CHKERRQ(ierr); 1317 v = a->a; 1318 for ( i=0; i<m; i++ ) { 1319 x = l[i]; 1320 M = a->i[i+1] - a->i[i]; 1321 for ( j=0; j<M; j++ ) { (*v++) *= x;} 1322 } 1323 ierr = VecRestoreArray(ll,&l); CHKERRQ(ierr); 1324 PLogFlops(nz); 1325 } 1326 if (rr) { 1327 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 1328 if (n != a->n) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Right scaling vector wrong length"); 1329 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1330 v = a->a; jj = a->j; 1331 for ( i=0; i<nz; i++ ) { 1332 (*v++) *= r[*jj++ + shift]; 1333 } 1334 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1335 PLogFlops(nz); 1336 } 1337 PetscFunctionReturn(0); 1338 } 1339 1340 #undef __FUNC__ 1341 #define __FUNC__ "MatGetSubMatrix_SeqAIJ" 1342 int MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,int csize,MatGetSubMatrixCall scall,Mat *B) 1343 { 1344 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data,*c; 1345 int *smap, i, k, kstart, kend, ierr, oldcols = a->n,*lens; 1346 int row,mat_i,*mat_j,tcol,first,step,*mat_ilen; 1347 register int sum,lensi; 1348 int *irow, *icol, nrows, ncols, shift = a->indexshift,*ssmap; 1349 int *starts,*j_new,*i_new,*aj = a->j, *ai = a->i,ii,*ailen = a->ilen; 1350 Scalar *a_new,*mat_a; 1351 Mat C; 1352 PetscTruth stride; 1353 1354 PetscFunctionBegin; 1355 ierr = ISSorted(isrow,(PetscTruth*)&i);CHKERRQ(ierr); 1356 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"ISrow is not sorted"); 1357 ierr = ISSorted(iscol,(PetscTruth*)&i);CHKERRQ(ierr); 1358 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"IScol is not sorted"); 1359 1360 ierr = ISGetIndices(isrow,&irow); CHKERRQ(ierr); 1361 ierr = ISGetSize(isrow,&nrows); CHKERRQ(ierr); 1362 ierr = ISGetSize(iscol,&ncols); CHKERRQ(ierr); 1363 1364 ierr = ISStrideGetInfo(iscol,&first,&step); CHKERRQ(ierr); 1365 ierr = ISStride(iscol,&stride); CHKERRQ(ierr); 1366 if (stride && step == 1) { 1367 /* special case of contiguous rows */ 1368 lens = (int *) PetscMalloc((ncols+nrows+1)*sizeof(int)); CHKPTRQ(lens); 1369 starts = lens + ncols; 1370 /* loop over new rows determining lens and starting points */ 1371 for (i=0; i<nrows; i++) { 1372 kstart = ai[irow[i]]+shift; 1373 kend = kstart + ailen[irow[i]]; 1374 for ( k=kstart; k<kend; k++ ) { 1375 if (aj[k]+shift >= first) { 1376 starts[i] = k; 1377 break; 1378 } 1379 } 1380 sum = 0; 1381 while (k < kend) { 1382 if (aj[k++]+shift >= first+ncols) break; 1383 sum++; 1384 } 1385 lens[i] = sum; 1386 } 1387 /* create submatrix */ 1388 if (scall == MAT_REUSE_MATRIX) { 1389 int n_cols,n_rows; 1390 ierr = MatGetSize(*B,&n_rows,&n_cols); CHKERRQ(ierr); 1391 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Reused submatrix wrong size"); 1392 ierr = MatZeroEntries(*B); CHKERRQ(ierr); 1393 C = *B; 1394 } else { 1395 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1396 } 1397 c = (Mat_SeqAIJ*) C->data; 1398 1399 /* loop over rows inserting into submatrix */ 1400 a_new = c->a; 1401 j_new = c->j; 1402 i_new = c->i; 1403 i_new[0] = -shift; 1404 for (i=0; i<nrows; i++) { 1405 ii = starts[i]; 1406 lensi = lens[i]; 1407 for ( k=0; k<lensi; k++ ) { 1408 *j_new++ = aj[ii+k] - first; 1409 } 1410 PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(Scalar)); 1411 a_new += lensi; 1412 i_new[i+1] = i_new[i] + lensi; 1413 c->ilen[i] = lensi; 1414 } 1415 PetscFree(lens); 1416 } else { 1417 ierr = ISGetIndices(iscol,&icol); CHKERRQ(ierr); 1418 smap = (int *) PetscMalloc((1+oldcols)*sizeof(int)); CHKPTRQ(smap); 1419 ssmap = smap + shift; 1420 lens = (int *) PetscMalloc((1+nrows)*sizeof(int)); CHKPTRQ(lens); 1421 PetscMemzero(smap,oldcols*sizeof(int)); 1422 for ( i=0; i<ncols; i++ ) smap[icol[i]] = i+1; 1423 /* determine lens of each row */ 1424 for (i=0; i<nrows; i++) { 1425 kstart = ai[irow[i]]+shift; 1426 kend = kstart + a->ilen[irow[i]]; 1427 lens[i] = 0; 1428 for ( k=kstart; k<kend; k++ ) { 1429 if (ssmap[aj[k]]) { 1430 lens[i]++; 1431 } 1432 } 1433 } 1434 /* Create and fill new matrix */ 1435 if (scall == MAT_REUSE_MATRIX) { 1436 c = (Mat_SeqAIJ *)((*B)->data); 1437 if (c->m != nrows || c->n != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Cannot reuse matrix. wrong size"); 1438 if (PetscMemcmp(c->ilen,lens, c->m *sizeof(int))) { 1439 SETERRQ(PETSC_ERR_ARG_SIZ,0,"Cannot reuse matrix. wrong no of nonzeros"); 1440 } 1441 PetscMemzero(c->ilen,c->m*sizeof(int)); 1442 C = *B; 1443 } else { 1444 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1445 } 1446 c = (Mat_SeqAIJ *)(C->data); 1447 for (i=0; i<nrows; i++) { 1448 row = irow[i]; 1449 kstart = ai[row]+shift; 1450 kend = kstart + a->ilen[row]; 1451 mat_i = c->i[i]+shift; 1452 mat_j = c->j + mat_i; 1453 mat_a = c->a + mat_i; 1454 mat_ilen = c->ilen + i; 1455 for ( k=kstart; k<kend; k++ ) { 1456 if ((tcol=ssmap[a->j[k]])) { 1457 *mat_j++ = tcol - (!shift); 1458 *mat_a++ = a->a[k]; 1459 (*mat_ilen)++; 1460 1461 } 1462 } 1463 } 1464 /* Free work space */ 1465 ierr = ISRestoreIndices(iscol,&icol); CHKERRQ(ierr); 1466 PetscFree(smap); PetscFree(lens); 1467 } 1468 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1469 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1470 1471 ierr = ISRestoreIndices(isrow,&irow); CHKERRQ(ierr); 1472 *B = C; 1473 PetscFunctionReturn(0); 1474 } 1475 1476 /* 1477 note: This can only work for identity for row and col. It would 1478 be good to check this and otherwise generate an error. 1479 */ 1480 #undef __FUNC__ 1481 #define __FUNC__ "MatILUFactor_SeqAIJ" 1482 int MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,double efill,int fill) 1483 { 1484 Mat_SeqAIJ *a = (Mat_SeqAIJ *) inA->data; 1485 int ierr; 1486 Mat outA; 1487 1488 PetscFunctionBegin; 1489 if (fill != 0) SETERRQ(PETSC_ERR_SUP,0,"Only fill=0 supported"); 1490 1491 outA = inA; 1492 inA->factor = FACTOR_LU; 1493 a->row = row; 1494 a->col = col; 1495 1496 /* Create the invert permutation so that it can be used in MatLUFactorNumeric() */ 1497 ierr = ISInvertPermutation(col,&(a->icol)); CHKERRQ(ierr); 1498 PLogObjectParent(inA,a->icol); 1499 1500 if (!a->solve_work) { /* this matrix may have been factored before */ 1501 a->solve_work = (Scalar *) PetscMalloc( (a->m+1)*sizeof(Scalar));CHKPTRQ(a->solve_work); 1502 } 1503 1504 if (!a->diag) { 1505 ierr = MatMarkDiag_SeqAIJ(inA); CHKERRQ(ierr); 1506 } 1507 ierr = MatLUFactorNumeric_SeqAIJ(inA,&outA); CHKERRQ(ierr); 1508 PetscFunctionReturn(0); 1509 } 1510 1511 #include "pinclude/blaslapack.h" 1512 #undef __FUNC__ 1513 #define __FUNC__ "MatScale_SeqAIJ" 1514 int MatScale_SeqAIJ(Scalar *alpha,Mat inA) 1515 { 1516 Mat_SeqAIJ *a = (Mat_SeqAIJ *) inA->data; 1517 int one = 1; 1518 1519 PetscFunctionBegin; 1520 BLscal_( &a->nz, alpha, a->a, &one ); 1521 PLogFlops(a->nz); 1522 PetscFunctionReturn(0); 1523 } 1524 1525 #undef __FUNC__ 1526 #define __FUNC__ "MatGetSubMatrices_SeqAIJ" 1527 int MatGetSubMatrices_SeqAIJ(Mat A,int n, IS *irow,IS *icol,MatGetSubMatrixCall scall,Mat **B) 1528 { 1529 int ierr,i; 1530 1531 PetscFunctionBegin; 1532 if (scall == MAT_INITIAL_MATRIX) { 1533 *B = (Mat *) PetscMalloc( (n+1)*sizeof(Mat) ); CHKPTRQ(*B); 1534 } 1535 1536 for ( i=0; i<n; i++ ) { 1537 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1538 } 1539 PetscFunctionReturn(0); 1540 } 1541 1542 #undef __FUNC__ 1543 #define __FUNC__ "MatGetBlockSize_SeqAIJ" 1544 int MatGetBlockSize_SeqAIJ(Mat A, int *bs) 1545 { 1546 PetscFunctionBegin; 1547 *bs = 1; 1548 PetscFunctionReturn(0); 1549 } 1550 1551 #undef __FUNC__ 1552 #define __FUNC__ "MatIncreaseOverlap_SeqAIJ" 1553 int MatIncreaseOverlap_SeqAIJ(Mat A, int is_max, IS *is, int ov) 1554 { 1555 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1556 int shift, row, i,j,k,l,m,n, *idx,ierr, *nidx, isz, val; 1557 int start, end, *ai, *aj; 1558 BT table; 1559 1560 PetscFunctionBegin; 1561 shift = a->indexshift; 1562 m = a->m; 1563 ai = a->i; 1564 aj = a->j+shift; 1565 1566 if (ov < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"illegal overlap value used"); 1567 1568 nidx = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(nidx); 1569 ierr = BTCreate(m,table); CHKERRQ(ierr); 1570 1571 for ( i=0; i<is_max; i++ ) { 1572 /* Initialize the two local arrays */ 1573 isz = 0; 1574 BTMemzero(m,table); 1575 1576 /* Extract the indices, assume there can be duplicate entries */ 1577 ierr = ISGetIndices(is[i],&idx); CHKERRQ(ierr); 1578 ierr = ISGetSize(is[i],&n); CHKERRQ(ierr); 1579 1580 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 1581 for ( j=0; j<n ; ++j){ 1582 if(!BTLookupSet(table, idx[j])) { nidx[isz++] = idx[j];} 1583 } 1584 ierr = ISRestoreIndices(is[i],&idx); CHKERRQ(ierr); 1585 ierr = ISDestroy(is[i]); CHKERRQ(ierr); 1586 1587 k = 0; 1588 for ( j=0; j<ov; j++){ /* for each overlap */ 1589 n = isz; 1590 for ( ; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 1591 row = nidx[k]; 1592 start = ai[row]; 1593 end = ai[row+1]; 1594 for ( l = start; l<end ; l++){ 1595 val = aj[l] + shift; 1596 if (!BTLookupSet(table,val)) {nidx[isz++] = val;} 1597 } 1598 } 1599 } 1600 ierr = ISCreateGeneral(PETSC_COMM_SELF, isz, nidx, (is+i)); CHKERRQ(ierr); 1601 } 1602 BTDestroy(table); 1603 PetscFree(nidx); 1604 PetscFunctionReturn(0); 1605 } 1606 1607 /* -------------------------------------------------------------- */ 1608 #undef __FUNC__ 1609 #define __FUNC__ "MatPermute_SeqAIJ" 1610 int MatPermute_SeqAIJ(Mat A, IS rowp, IS colp, Mat *B) 1611 { 1612 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1613 Scalar *vwork; 1614 int i, ierr, nz, m = a->m, n = a->n, *cwork; 1615 int *row,*col,*cnew,j,*lens; 1616 IS icolp,irowp; 1617 1618 PetscFunctionBegin; 1619 ierr = ISInvertPermutation(rowp,&irowp); CHKERRQ(ierr); 1620 ierr = ISGetIndices(irowp,&row); CHKERRQ(ierr); 1621 ierr = ISInvertPermutation(colp,&icolp); CHKERRQ(ierr); 1622 ierr = ISGetIndices(icolp,&col); CHKERRQ(ierr); 1623 1624 /* determine lengths of permuted rows */ 1625 lens = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(lens); 1626 for (i=0; i<m; i++ ) { 1627 lens[row[i]] = a->i[i+1] - a->i[i]; 1628 } 1629 ierr = MatCreateSeqAIJ(A->comm,m,n,0,lens,B);CHKERRQ(ierr); 1630 PetscFree(lens); 1631 1632 cnew = (int *) PetscMalloc( n*sizeof(int) ); CHKPTRQ(cnew); 1633 for (i=0; i<m; i++) { 1634 ierr = MatGetRow(A,i,&nz,&cwork,&vwork); CHKERRQ(ierr); 1635 for (j=0; j<nz; j++ ) { cnew[j] = col[cwork[j]];} 1636 ierr = MatSetValues(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES); CHKERRQ(ierr); 1637 ierr = MatRestoreRow(A,i,&nz,&cwork,&vwork); CHKERRQ(ierr); 1638 } 1639 PetscFree(cnew); 1640 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1641 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1642 ierr = ISRestoreIndices(irowp,&row); CHKERRQ(ierr); 1643 ierr = ISRestoreIndices(icolp,&col); CHKERRQ(ierr); 1644 ierr = ISDestroy(irowp); CHKERRQ(ierr); 1645 ierr = ISDestroy(icolp); CHKERRQ(ierr); 1646 PetscFunctionReturn(0); 1647 } 1648 1649 #undef __FUNC__ 1650 #define __FUNC__ "MatPrintHelp_SeqAIJ" 1651 int MatPrintHelp_SeqAIJ(Mat A) 1652 { 1653 static int called = 0; 1654 MPI_Comm comm = A->comm; 1655 1656 PetscFunctionBegin; 1657 if (called) {PetscFunctionReturn(0);} else called = 1; 1658 (*PetscHelpPrintf)(comm," Options for MATSEQAIJ and MATMPIAIJ matrix formats (the defaults):\n"); 1659 (*PetscHelpPrintf)(comm," -mat_lu_pivotthreshold <threshold>: Set pivoting threshold\n"); 1660 (*PetscHelpPrintf)(comm," -mat_aij_oneindex: internal indices begin at 1 instead of the default 0.\n"); 1661 (*PetscHelpPrintf)(comm," -mat_aij_no_inode: Do not use inodes\n"); 1662 (*PetscHelpPrintf)(comm," -mat_aij_inode_limit <limit>: Set inode limit (max limit=5)\n"); 1663 #if defined(HAVE_ESSL) 1664 (*PetscHelpPrintf)(comm," -mat_aij_essl: Use IBM sparse LU factorization and solve.\n"); 1665 #endif 1666 PetscFunctionReturn(0); 1667 } 1668 extern int MatEqual_SeqAIJ(Mat A,Mat B, PetscTruth* flg); 1669 extern int MatFDColoringCreate_SeqAIJ(Mat,ISColoring,MatFDColoring); 1670 extern int MatColoringPatch_SeqAIJ(Mat,int,int *,ISColoring *); 1671 1672 /* -------------------------------------------------------------------*/ 1673 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 1674 MatGetRow_SeqAIJ,MatRestoreRow_SeqAIJ, 1675 MatMult_SeqAIJ,MatMultAdd_SeqAIJ, 1676 MatMultTrans_SeqAIJ,MatMultTransAdd_SeqAIJ, 1677 MatSolve_SeqAIJ,MatSolveAdd_SeqAIJ, 1678 MatSolveTrans_SeqAIJ,MatSolveTransAdd_SeqAIJ, 1679 MatLUFactor_SeqAIJ,0, 1680 MatRelax_SeqAIJ, 1681 MatTranspose_SeqAIJ, 1682 MatGetInfo_SeqAIJ,MatEqual_SeqAIJ, 1683 MatGetDiagonal_SeqAIJ,MatDiagonalScale_SeqAIJ,MatNorm_SeqAIJ, 1684 0,MatAssemblyEnd_SeqAIJ, 1685 MatCompress_SeqAIJ, 1686 MatSetOption_SeqAIJ,MatZeroEntries_SeqAIJ,MatZeroRows_SeqAIJ, 1687 MatLUFactorSymbolic_SeqAIJ,MatLUFactorNumeric_SeqAIJ,0,0, 1688 MatGetSize_SeqAIJ,MatGetSize_SeqAIJ,MatGetOwnershipRange_SeqAIJ, 1689 MatILUFactorSymbolic_SeqAIJ,0, 1690 0,0, 1691 MatConvertSameType_SeqAIJ,0,0, 1692 MatILUFactor_SeqAIJ,0,0, 1693 MatGetSubMatrices_SeqAIJ,MatIncreaseOverlap_SeqAIJ, 1694 MatGetValues_SeqAIJ,0, 1695 MatPrintHelp_SeqAIJ, 1696 MatScale_SeqAIJ,0,0, 1697 MatILUDTFactor_SeqAIJ, 1698 MatGetBlockSize_SeqAIJ, 1699 MatGetRowIJ_SeqAIJ, 1700 MatRestoreRowIJ_SeqAIJ, 1701 MatGetColumnIJ_SeqAIJ, 1702 MatRestoreColumnIJ_SeqAIJ, 1703 MatFDColoringCreate_SeqAIJ, 1704 MatColoringPatch_SeqAIJ, 1705 0, 1706 MatPermute_SeqAIJ, 1707 0, 1708 0, 1709 0, 1710 0, 1711 MatGetMaps_Petsc}; 1712 1713 extern int MatUseSuperLU_SeqAIJ(Mat); 1714 extern int MatUseEssl_SeqAIJ(Mat); 1715 extern int MatUseDXML_SeqAIJ(Mat); 1716 1717 #undef __FUNC__ 1718 #define __FUNC__ "MatSeqAIJSetColumnIndices_SeqAIJ" 1719 int MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,int *indices) 1720 { 1721 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 1722 int i,nz,n; 1723 1724 PetscFunctionBegin; 1725 if (aij->indexshift) SETERRQ(1,1,"No support with 1 based indexing"); 1726 1727 nz = aij->maxnz; 1728 n = aij->n; 1729 for (i=0; i<nz; i++) { 1730 aij->j[i] = indices[i]; 1731 } 1732 aij->nz = nz; 1733 for ( i=0; i<n; i++ ) { 1734 aij->ilen[i] = aij->imax[i]; 1735 } 1736 1737 PetscFunctionReturn(0); 1738 } 1739 1740 #undef __FUNC__ 1741 #define __FUNC__ "MatSeqAIJSetColumnIndices" 1742 /*@ 1743 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 1744 in the matrix. 1745 1746 Input Parameters: 1747 + mat - the SeqAIJ matrix 1748 - indices - the column indices 1749 1750 Notes: 1751 This can be called if you have precomputed the nonzero structure of the 1752 matrix and want to provide it to the matrix object to improve the performance 1753 of the MatSetValues() operation. 1754 1755 You MUST have set the correct numbers of nonzeros per row in the call to 1756 MatCreateSeqAIJ(). 1757 1758 MUST be called before any calls to MatSetValues(); 1759 1760 @*/ 1761 int MatSeqAIJSetColumnIndices(Mat mat,int *indices) 1762 { 1763 int ierr,(*f)(Mat,int *); 1764 1765 PetscFunctionBegin; 1766 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1767 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqAIJSetColumnIndices_C",(void **)&f); 1768 CHKERRQ(ierr); 1769 if (f) { 1770 ierr = (*f)(mat,indices);CHKERRQ(ierr); 1771 } else { 1772 SETERRQ(1,1,"Wrong type of matrix to set column indices"); 1773 } 1774 PetscFunctionReturn(0); 1775 } 1776 1777 #undef __FUNC__ 1778 #define __FUNC__ "MatCreateSeqAIJ" 1779 /*@C 1780 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 1781 (the default parallel PETSc format). For good matrix assembly performance 1782 the user should preallocate the matrix storage by setting the parameter nz 1783 (or the array nzz). By setting these parameters accurately, performance 1784 during matrix assembly can be increased by more than a factor of 50. 1785 1786 Collective on MPI_Comm 1787 1788 Input Parameters: 1789 + comm - MPI communicator, set to PETSC_COMM_SELF 1790 . m - number of rows 1791 . n - number of columns 1792 . nz - number of nonzeros per row (same for all rows) 1793 - nzz - array containing the number of nonzeros in the various rows 1794 (possibly different for each row) or PETSC_NULL 1795 1796 Output Parameter: 1797 . A - the matrix 1798 1799 Notes: 1800 The AIJ format (also called the Yale sparse matrix format or 1801 compressed row storage), is fully compatible with standard Fortran 77 1802 storage. That is, the stored row and column indices can begin at 1803 either one (as in Fortran) or zero. See the users' manual for details. 1804 1805 Specify the preallocated storage with either nz or nnz (not both). 1806 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1807 allocation. For large problems you MUST preallocate memory or you 1808 will get TERRIBLE performance, see the users' manual chapter on matrices. 1809 1810 By default, this format uses inodes (identical nodes) when possible, to 1811 improve numerical efficiency of matrix-vector products and solves. We 1812 search for consecutive rows with the same nonzero structure, thereby 1813 reusing matrix information to achieve increased efficiency. 1814 1815 Options Database Keys: 1816 + -mat_aij_no_inode - Do not use inodes 1817 . -mat_aij_inode_limit <limit> - Sets inode limit (max limit=5) 1818 - -mat_aij_oneindex - Internally use indexing starting at 1 1819 rather than 0. Note that when calling MatSetValues(), 1820 the user still MUST index entries starting at 0! 1821 1822 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices() 1823 @*/ 1824 int MatCreateSeqAIJ(MPI_Comm comm,int m,int n,int nz,int *nnz, Mat *A) 1825 { 1826 Mat B; 1827 Mat_SeqAIJ *b; 1828 int i, len, ierr, flg,size; 1829 1830 PetscFunctionBegin; 1831 MPI_Comm_size(comm,&size); 1832 if (size > 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Comm must be of size 1"); 1833 1834 *A = 0; 1835 PetscHeaderCreate(B,_p_Mat,struct _MatOps,MAT_COOKIE,MATSEQAIJ,comm,MatDestroy,MatView); 1836 PLogObjectCreate(B); 1837 B->data = (void *) (b = PetscNew(Mat_SeqAIJ)); CHKPTRQ(b); 1838 PetscMemzero(b,sizeof(Mat_SeqAIJ)); 1839 PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps)); 1840 B->ops->destroy = MatDestroy_SeqAIJ; 1841 B->ops->view = MatView_SeqAIJ; 1842 B->factor = 0; 1843 B->lupivotthreshold = 1.0; 1844 B->mapping = 0; 1845 ierr = OptionsGetDouble(PETSC_NULL,"-mat_lu_pivotthreshold",&B->lupivotthreshold,&flg); CHKERRQ(ierr); 1846 b->ilu_preserve_row_sums = PETSC_FALSE; 1847 ierr = OptionsHasName(PETSC_NULL,"-pc_ilu_preserve_row_sums",(int*) &b->ilu_preserve_row_sums);CHKERRQ(ierr); 1848 b->row = 0; 1849 b->col = 0; 1850 b->icol = 0; 1851 b->indexshift = 0; 1852 b->reallocs = 0; 1853 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_oneindex", &flg); CHKERRQ(ierr); 1854 if (flg) b->indexshift = -1; 1855 1856 b->m = m; B->m = m; B->M = m; 1857 b->n = n; B->n = n; B->N = n; 1858 b->imax = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(b->imax); 1859 if (nnz == PETSC_NULL) { 1860 if (nz == PETSC_DEFAULT) nz = 10; 1861 else if (nz <= 0) nz = 1; 1862 for ( i=0; i<m; i++ ) b->imax[i] = nz; 1863 nz = nz*m; 1864 } else { 1865 nz = 0; 1866 for ( i=0; i<m; i++ ) {b->imax[i] = nnz[i]; nz += nnz[i];} 1867 } 1868 1869 /* allocate the matrix space */ 1870 len = nz*(sizeof(int) + sizeof(Scalar)) + (b->m+1)*sizeof(int); 1871 b->a = (Scalar *) PetscMalloc( len ); CHKPTRQ(b->a); 1872 b->j = (int *) (b->a + nz); 1873 PetscMemzero(b->j,nz*sizeof(int)); 1874 b->i = b->j + nz; 1875 b->singlemalloc = 1; 1876 1877 b->i[0] = -b->indexshift; 1878 for (i=1; i<m+1; i++) { 1879 b->i[i] = b->i[i-1] + b->imax[i-1]; 1880 } 1881 1882 /* b->ilen will count nonzeros in each row so far. */ 1883 b->ilen = (int *) PetscMalloc((m+1)*sizeof(int)); 1884 PLogObjectMemory(B,len+2*(m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 1885 for ( i=0; i<b->m; i++ ) { b->ilen[i] = 0;} 1886 1887 b->nz = 0; 1888 b->maxnz = nz; 1889 b->sorted = 0; 1890 b->roworiented = 1; 1891 b->nonew = 0; 1892 b->diag = 0; 1893 b->solve_work = 0; 1894 b->spptr = 0; 1895 b->inode.node_count = 0; 1896 b->inode.size = 0; 1897 b->inode.limit = 5; 1898 b->inode.max_limit = 5; 1899 B->info.nz_unneeded = (double)b->maxnz; 1900 1901 *A = B; 1902 1903 /* SuperLU is not currently supported through PETSc */ 1904 #if defined(HAVE_SUPERLU) 1905 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_superlu", &flg); CHKERRQ(ierr); 1906 if (flg) { ierr = MatUseSuperLU_SeqAIJ(B); CHKERRQ(ierr); } 1907 #endif 1908 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_essl", &flg); CHKERRQ(ierr); 1909 if (flg) { ierr = MatUseEssl_SeqAIJ(B); CHKERRQ(ierr); } 1910 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_dxml", &flg); CHKERRQ(ierr); 1911 if (flg) { 1912 if (!b->indexshift) SETERRQ( PETSC_ERR_LIB,0,"need -mat_aij_oneindex with -mat_aij_dxml"); 1913 ierr = MatUseDXML_SeqAIJ(B); CHKERRQ(ierr); 1914 } 1915 ierr = OptionsHasName(PETSC_NULL,"-help", &flg); CHKERRQ(ierr); 1916 if (flg) {ierr = MatPrintHelp(B); CHKERRQ(ierr); } 1917 1918 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 1919 "MatSeqAIJSetColumnIndices_SeqAIJ", 1920 (void*)MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 1921 PetscFunctionReturn(0); 1922 } 1923 1924 #undef __FUNC__ 1925 #define __FUNC__ "MatConvertSameType_SeqAIJ" 1926 int MatConvertSameType_SeqAIJ(Mat A,Mat *B,int cpvalues) 1927 { 1928 Mat C; 1929 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ *) A->data; 1930 int i,len, m = a->m,shift = a->indexshift,ierr; 1931 1932 PetscFunctionBegin; 1933 *B = 0; 1934 PetscHeaderCreate(C,_p_Mat,struct _MatOps,MAT_COOKIE,MATSEQAIJ,A->comm,MatDestroy,MatView); 1935 PLogObjectCreate(C); 1936 C->data = (void *) (c = PetscNew(Mat_SeqAIJ)); CHKPTRQ(c); 1937 PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps)); 1938 C->ops->destroy = MatDestroy_SeqAIJ; 1939 C->ops->view = MatView_SeqAIJ; 1940 C->factor = A->factor; 1941 c->row = 0; 1942 c->col = 0; 1943 c->icol = 0; 1944 c->indexshift = shift; 1945 C->assembled = PETSC_TRUE; 1946 1947 c->m = C->m = a->m; 1948 c->n = C->n = a->n; 1949 C->M = a->m; 1950 C->N = a->n; 1951 1952 c->imax = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(c->imax); 1953 c->ilen = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(c->ilen); 1954 for ( i=0; i<m; i++ ) { 1955 c->imax[i] = a->imax[i]; 1956 c->ilen[i] = a->ilen[i]; 1957 } 1958 1959 /* allocate the matrix space */ 1960 c->singlemalloc = 1; 1961 len = (m+1)*sizeof(int)+(a->i[m])*(sizeof(Scalar)+sizeof(int)); 1962 c->a = (Scalar *) PetscMalloc( len ); CHKPTRQ(c->a); 1963 c->j = (int *) (c->a + a->i[m] + shift); 1964 c->i = c->j + a->i[m] + shift; 1965 PetscMemcpy(c->i,a->i,(m+1)*sizeof(int)); 1966 if (m > 0) { 1967 PetscMemcpy(c->j,a->j,(a->i[m]+shift)*sizeof(int)); 1968 if (cpvalues == COPY_VALUES) { 1969 PetscMemcpy(c->a,a->a,(a->i[m]+shift)*sizeof(Scalar)); 1970 } 1971 } 1972 1973 PLogObjectMemory(C,len+2*(m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 1974 c->sorted = a->sorted; 1975 c->roworiented = a->roworiented; 1976 c->nonew = a->nonew; 1977 c->ilu_preserve_row_sums = a->ilu_preserve_row_sums; 1978 1979 if (a->diag) { 1980 c->diag = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(c->diag); 1981 PLogObjectMemory(C,(m+1)*sizeof(int)); 1982 for ( i=0; i<m; i++ ) { 1983 c->diag[i] = a->diag[i]; 1984 } 1985 } else c->diag = 0; 1986 c->inode.limit = a->inode.limit; 1987 c->inode.max_limit = a->inode.max_limit; 1988 if (a->inode.size){ 1989 c->inode.size = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(c->inode.size); 1990 c->inode.node_count = a->inode.node_count; 1991 PetscMemcpy( c->inode.size, a->inode.size, (m+1)*sizeof(int)); 1992 } else { 1993 c->inode.size = 0; 1994 c->inode.node_count = 0; 1995 } 1996 c->nz = a->nz; 1997 c->maxnz = a->maxnz; 1998 c->solve_work = 0; 1999 c->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 2000 2001 *B = C; 2002 ierr = PetscObjectComposeFunction((PetscObject)C,"MatSeqAIJSetColumnIndices_C", 2003 "MatSeqAIJSetColumnIndices_SeqAIJ", 2004 (void*)MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 2005 PetscFunctionReturn(0); 2006 } 2007 2008 #undef __FUNC__ 2009 #define __FUNC__ "MatLoad_SeqAIJ" 2010 int MatLoad_SeqAIJ(Viewer viewer,MatType type,Mat *A) 2011 { 2012 Mat_SeqAIJ *a; 2013 Mat B; 2014 int i, nz, ierr, fd, header[4],size,*rowlengths = 0,M,N,shift; 2015 MPI_Comm comm; 2016 2017 PetscFunctionBegin; 2018 PetscObjectGetComm((PetscObject) viewer,&comm); 2019 MPI_Comm_size(comm,&size); 2020 if (size > 1) SETERRQ(PETSC_ERR_ARG_SIZ,0,"view must have one processor"); 2021 ierr = ViewerBinaryGetDescriptor(viewer,&fd); CHKERRQ(ierr); 2022 ierr = PetscBinaryRead(fd,header,4,PETSC_INT); CHKERRQ(ierr); 2023 if (header[0] != MAT_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,0,"not matrix object in file"); 2024 M = header[1]; N = header[2]; nz = header[3]; 2025 2026 if (nz < 0) { 2027 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,1,"Matrix stored in special format on disk, cannot load as SeqAIJ"); 2028 } 2029 2030 /* read in row lengths */ 2031 rowlengths = (int*) PetscMalloc( M*sizeof(int) ); CHKPTRQ(rowlengths); 2032 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT); CHKERRQ(ierr); 2033 2034 /* create our matrix */ 2035 ierr = MatCreateSeqAIJ(comm,M,N,0,rowlengths,A); CHKERRQ(ierr); 2036 B = *A; 2037 a = (Mat_SeqAIJ *) B->data; 2038 shift = a->indexshift; 2039 2040 /* read in column indices and adjust for Fortran indexing*/ 2041 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT); CHKERRQ(ierr); 2042 if (shift) { 2043 for ( i=0; i<nz; i++ ) { 2044 a->j[i] += 1; 2045 } 2046 } 2047 2048 /* read in nonzero values */ 2049 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR); CHKERRQ(ierr); 2050 2051 /* set matrix "i" values */ 2052 a->i[0] = -shift; 2053 for ( i=1; i<= M; i++ ) { 2054 a->i[i] = a->i[i-1] + rowlengths[i-1]; 2055 a->ilen[i-1] = rowlengths[i-1]; 2056 } 2057 PetscFree(rowlengths); 2058 2059 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 2060 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 2061 PetscFunctionReturn(0); 2062 } 2063 2064 #undef __FUNC__ 2065 #define __FUNC__ "MatEqual_SeqAIJ" 2066 int MatEqual_SeqAIJ(Mat A,Mat B, PetscTruth* flg) 2067 { 2068 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)B->data; 2069 2070 PetscFunctionBegin; 2071 if (B->type !=MATSEQAIJ)SETERRQ(PETSC_ERR_ARG_INCOMP,0,"Matrices must be same type"); 2072 2073 /* If the matrix dimensions are not equal, or no of nonzeros or shift */ 2074 if ((a->m != b->m ) || (a->n !=b->n) ||( a->nz != b->nz)|| 2075 (a->indexshift != b->indexshift)) { 2076 *flg = PETSC_FALSE; PetscFunctionReturn(0); 2077 } 2078 2079 /* if the a->i are the same */ 2080 if (PetscMemcmp(a->i,b->i,(a->m+1)*sizeof(int))) { 2081 *flg = PETSC_FALSE; PetscFunctionReturn(0); 2082 } 2083 2084 /* if a->j are the same */ 2085 if (PetscMemcmp(a->j, b->j, (a->nz)*sizeof(int))) { 2086 *flg = PETSC_FALSE; PetscFunctionReturn(0); 2087 } 2088 2089 /* if a->a are the same */ 2090 if (PetscMemcmp(a->a, b->a, (a->nz)*sizeof(Scalar))) { 2091 *flg = PETSC_FALSE; PetscFunctionReturn(0); 2092 } 2093 *flg = PETSC_TRUE; 2094 PetscFunctionReturn(0); 2095 2096 } 2097