1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: aij.c,v 1.282 1998/09/25 03:14:31 bsmith 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 "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 (A->rmap) { 671 ierr = MapDestroy(A->rmap);CHKERRQ(ierr); 672 } 673 if (A->cmap) { 674 ierr = MapDestroy(A->cmap);CHKERRQ(ierr); 675 } 676 #if defined(USE_PETSC_LOG) 677 PLogObjectState((PetscObject)A,"Rows=%d, Cols=%d, NZ=%d",a->m,a->n,a->nz); 678 #endif 679 PetscFree(a->a); 680 if (!a->singlemalloc) { PetscFree(a->i); PetscFree(a->j);} 681 if (a->diag) PetscFree(a->diag); 682 if (a->ilen) PetscFree(a->ilen); 683 if (a->imax) PetscFree(a->imax); 684 if (a->solve_work) PetscFree(a->solve_work); 685 if (a->inode.size) PetscFree(a->inode.size); 686 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} 687 if (a->saved_values) PetscFree(a->saved_values); 688 PetscFree(a); 689 690 PLogObjectDestroy(A); 691 PetscHeaderDestroy(A); 692 PetscFunctionReturn(0); 693 } 694 695 #undef __FUNC__ 696 #define __FUNC__ "MatCompress_SeqAIJ" 697 int MatCompress_SeqAIJ(Mat A) 698 { 699 PetscFunctionBegin; 700 PetscFunctionReturn(0); 701 } 702 703 #undef __FUNC__ 704 #define __FUNC__ "MatSetOption_SeqAIJ" 705 int MatSetOption_SeqAIJ(Mat A,MatOption op) 706 { 707 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 708 709 PetscFunctionBegin; 710 if (op == MAT_ROW_ORIENTED) a->roworiented = 1; 711 else if (op == MAT_COLUMN_ORIENTED) a->roworiented = 0; 712 else if (op == MAT_COLUMNS_SORTED) a->sorted = 1; 713 else if (op == MAT_COLUMNS_UNSORTED) a->sorted = 0; 714 else if (op == MAT_NO_NEW_NONZERO_LOCATIONS) a->nonew = 1; 715 else if (op == MAT_NEW_NONZERO_LOCATION_ERROR) a->nonew = -1; 716 else if (op == MAT_NEW_NONZERO_ALLOCATION_ERROR) a->nonew = -2; 717 else if (op == MAT_YES_NEW_NONZERO_LOCATIONS) a->nonew = 0; 718 else if (op == MAT_ROWS_SORTED || 719 op == MAT_ROWS_UNSORTED || 720 op == MAT_SYMMETRIC || 721 op == MAT_STRUCTURALLY_SYMMETRIC || 722 op == MAT_YES_NEW_DIAGONALS || 723 op == MAT_IGNORE_OFF_PROC_ENTRIES|| 724 op == MAT_USE_HASH_TABLE) 725 PLogInfo(A,"MatSetOption_SeqAIJ:Option ignored\n"); 726 else if (op == MAT_NO_NEW_DIAGONALS) { 727 SETERRQ(PETSC_ERR_SUP,0,"MAT_NO_NEW_DIAGONALS"); 728 } else if (op == MAT_INODE_LIMIT_1) a->inode.limit = 1; 729 else if (op == MAT_INODE_LIMIT_2) a->inode.limit = 2; 730 else if (op == MAT_INODE_LIMIT_3) a->inode.limit = 3; 731 else if (op == MAT_INODE_LIMIT_4) a->inode.limit = 4; 732 else if (op == MAT_INODE_LIMIT_5) a->inode.limit = 5; 733 else SETERRQ(PETSC_ERR_SUP,0,"unknown option"); 734 PetscFunctionReturn(0); 735 } 736 737 #undef __FUNC__ 738 #define __FUNC__ "MatGetDiagonal_SeqAIJ" 739 int MatGetDiagonal_SeqAIJ(Mat A,Vec v) 740 { 741 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 742 int i,j, n,shift = a->indexshift,ierr; 743 Scalar *x, zero = 0.0; 744 745 PetscFunctionBegin; 746 ierr = VecSet(&zero,v);CHKERRQ(ierr); 747 ierr = VecGetArray(v,&x);;CHKERRQ(ierr); 748 ierr = VecGetLocalSize(v,&n);;CHKERRQ(ierr); 749 if (n != a->m) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Nonconforming matrix and vector"); 750 for ( i=0; i<a->m; i++ ) { 751 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 752 if (a->j[j]+shift == i) { 753 x[i] = a->a[j]; 754 break; 755 } 756 } 757 } 758 ierr = VecRestoreArray(v,&x);;CHKERRQ(ierr); 759 PetscFunctionReturn(0); 760 } 761 762 /* -------------------------------------------------------*/ 763 /* Should check that shapes of vectors and matrices match */ 764 /* -------------------------------------------------------*/ 765 #undef __FUNC__ 766 #define __FUNC__ "MatMultTrans_SeqAIJ" 767 int MatMultTrans_SeqAIJ(Mat A,Vec xx,Vec yy) 768 { 769 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 770 Scalar *x, *y, *v, alpha; 771 int ierr,m = a->m, n, i, *idx, shift = a->indexshift; 772 773 PetscFunctionBegin; 774 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 775 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 776 PetscMemzero(y,a->n*sizeof(Scalar)); 777 y = y + shift; /* shift for Fortran start by 1 indexing */ 778 for ( i=0; i<m; i++ ) { 779 idx = a->j + a->i[i] + shift; 780 v = a->a + a->i[i] + shift; 781 n = a->i[i+1] - a->i[i]; 782 alpha = x[i]; 783 while (n-->0) {y[*idx++] += alpha * *v++;} 784 } 785 PLogFlops(2*a->nz - a->n); 786 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 787 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 788 PetscFunctionReturn(0); 789 } 790 791 #undef __FUNC__ 792 #define __FUNC__ "MatMultTransAdd_SeqAIJ" 793 int MatMultTransAdd_SeqAIJ(Mat A,Vec xx,Vec zz,Vec yy) 794 { 795 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 796 Scalar *x, *y, *v, alpha; 797 int ierr,m = a->m, n, i, *idx,shift = a->indexshift; 798 799 PetscFunctionBegin; 800 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 801 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 802 if (zz != yy) VecCopy(zz,yy); 803 y = y + shift; /* shift for Fortran start by 1 indexing */ 804 for ( i=0; i<m; i++ ) { 805 idx = a->j + a->i[i] + shift; 806 v = a->a + a->i[i] + shift; 807 n = a->i[i+1] - a->i[i]; 808 alpha = x[i]; 809 while (n-->0) {y[*idx++] += alpha * *v++;} 810 } 811 PLogFlops(2*a->nz); 812 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 813 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 814 PetscFunctionReturn(0); 815 } 816 817 #undef __FUNC__ 818 #define __FUNC__ "MatMult_SeqAIJ" 819 int MatMult_SeqAIJ(Mat A,Vec xx,Vec yy) 820 { 821 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 822 Scalar *x, *y, *v, sum; 823 int ierr,m = a->m, *idx, shift = a->indexshift,*ii; 824 #if !defined(USE_FORTRAN_KERNEL_MULTAIJ) 825 int n, i, jrow,j; 826 #endif 827 828 #if defined(HAVE_PRAGMA_DISJOINT) 829 #pragma disjoint(*x,*y,*v) 830 #endif 831 832 PetscFunctionBegin; 833 ierr = VecGetArray(xx,&x); CHKERRQ(ierr); 834 ierr = VecGetArray(yy,&y); CHKERRQ(ierr); 835 x = x + shift; /* shift for Fortran start by 1 indexing */ 836 idx = a->j; 837 v = a->a; 838 ii = a->i; 839 #if defined(USE_FORTRAN_KERNEL_MULTAIJ) 840 fortranmultaij_(&m,x,ii,idx+shift,v+shift,y); 841 #else 842 v += shift; /* shift for Fortran start by 1 indexing */ 843 idx += shift; 844 for ( i=0; i<m; i++ ) { 845 jrow = ii[i]; 846 n = ii[i+1] - jrow; 847 sum = 0.0; 848 for ( j=0; j<n; j++) { 849 sum += v[jrow]*x[idx[jrow]]; jrow++; 850 } 851 y[i] = sum; 852 } 853 #endif 854 PLogFlops(2*a->nz - m); 855 ierr = VecRestoreArray(xx,&x); CHKERRQ(ierr); 856 ierr = VecRestoreArray(yy,&y); CHKERRQ(ierr); 857 PetscFunctionReturn(0); 858 } 859 860 #undef __FUNC__ 861 #define __FUNC__ "MatMultAdd_SeqAIJ" 862 int MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 863 { 864 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 865 Scalar *x, *y, *z, *v, sum; 866 int ierr,m = a->m, *idx, shift = a->indexshift,*ii; 867 #if !defined(USE_FORTRAN_KERNEL_MULTADDAIJ) 868 int n,i,jrow,j; 869 #endif 870 871 PetscFunctionBegin; 872 ierr = VecGetArray(xx,&x); CHKERRQ(ierr); 873 ierr = VecGetArray(yy,&y); CHKERRQ(ierr); 874 ierr = VecGetArray(zz,&z); CHKERRQ(ierr); 875 x = x + shift; /* shift for Fortran start by 1 indexing */ 876 idx = a->j; 877 v = a->a; 878 ii = a->i; 879 #if defined(USE_FORTRAN_KERNEL_MULTADDAIJ) 880 fortranmultaddaij_(&m,x,ii,idx+shift,v+shift,y,z); 881 #else 882 v += shift; /* shift for Fortran start by 1 indexing */ 883 idx += shift; 884 for ( i=0; i<m; i++ ) { 885 jrow = ii[i]; 886 n = ii[i+1] - jrow; 887 sum = y[i]; 888 for ( j=0; j<n; j++) { 889 sum += v[jrow]*x[idx[jrow]]; jrow++; 890 } 891 z[i] = sum; 892 } 893 #endif 894 PLogFlops(2*a->nz); 895 ierr = VecRestoreArray(xx,&x); CHKERRQ(ierr); 896 ierr = VecRestoreArray(yy,&y); CHKERRQ(ierr); 897 ierr = VecRestoreArray(zz,&z); CHKERRQ(ierr); 898 PetscFunctionReturn(0); 899 } 900 901 /* 902 Adds diagonal pointers to sparse matrix structure. 903 */ 904 905 #undef __FUNC__ 906 #define __FUNC__ "MatMarkDiag_SeqAIJ" 907 int MatMarkDiag_SeqAIJ(Mat A) 908 { 909 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 910 int i,j, *diag, m = a->m,shift = a->indexshift; 911 912 PetscFunctionBegin; 913 diag = (int *) PetscMalloc( (m+1)*sizeof(int)); CHKPTRQ(diag); 914 PLogObjectMemory(A,(m+1)*sizeof(int)); 915 for ( i=0; i<a->m; i++ ) { 916 diag[i] = a->i[i+1]; 917 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 918 if (a->j[j]+shift == i) { 919 diag[i] = j - shift; 920 break; 921 } 922 } 923 } 924 a->diag = diag; 925 PetscFunctionReturn(0); 926 } 927 928 #undef __FUNC__ 929 #define __FUNC__ "MatRelax_SeqAIJ" 930 int MatRelax_SeqAIJ(Mat A,Vec bb,double omega,MatSORType flag, 931 double fshift,int its,Vec xx) 932 { 933 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 934 Scalar *x, *b, *bs, d, *xs, sum, *v = a->a,*t,scale,*ts, *xb; 935 int ierr, *idx, *diag,n = a->n, m = a->m, i, shift = a->indexshift; 936 937 PetscFunctionBegin; 938 ierr = VecGetArray(xx,&x); CHKERRQ(ierr); 939 ierr = VecGetArray(bb,&b);CHKERRQ(ierr); 940 if (!a->diag) {ierr = MatMarkDiag_SeqAIJ(A);CHKERRQ(ierr);} 941 diag = a->diag; 942 xs = x + shift; /* shifted by one for index start of a or a->j*/ 943 if (flag == SOR_APPLY_UPPER) { 944 /* apply ( U + D/omega) to the vector */ 945 bs = b + shift; 946 for ( i=0; i<m; i++ ) { 947 d = fshift + a->a[diag[i] + shift]; 948 n = a->i[i+1] - diag[i] - 1; 949 PLogFlops(2*n-1); 950 idx = a->j + diag[i] + (!shift); 951 v = a->a + diag[i] + (!shift); 952 sum = b[i]*d/omega; 953 SPARSEDENSEDOT(sum,bs,v,idx,n); 954 x[i] = sum; 955 } 956 ierr = VecRestoreArray(xx,&x); CHKERRQ(ierr); 957 ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr); 958 PetscFunctionReturn(0); 959 } 960 if (flag == SOR_APPLY_LOWER) { 961 SETERRQ(PETSC_ERR_SUP,0,"SOR_APPLY_LOWER is not done"); 962 } else if (flag & SOR_EISENSTAT) { 963 /* Let A = L + U + D; where L is lower trianglar, 964 U is upper triangular, E is diagonal; This routine applies 965 966 (L + E)^{-1} A (U + E)^{-1} 967 968 to a vector efficiently using Eisenstat's trick. This is for 969 the case of SSOR preconditioner, so E is D/omega where omega 970 is the relaxation factor. 971 */ 972 t = (Scalar *) PetscMalloc( m*sizeof(Scalar) ); CHKPTRQ(t); 973 scale = (2.0/omega) - 1.0; 974 975 /* x = (E + U)^{-1} b */ 976 for ( i=m-1; i>=0; i-- ) { 977 d = fshift + a->a[diag[i] + shift]; 978 n = a->i[i+1] - diag[i] - 1; 979 PLogFlops(2*n-1); 980 idx = a->j + diag[i] + (!shift); 981 v = a->a + diag[i] + (!shift); 982 sum = b[i]; 983 SPARSEDENSEMDOT(sum,xs,v,idx,n); 984 x[i] = omega*(sum/d); 985 } 986 987 /* t = b - (2*E - D)x */ 988 v = a->a; 989 for ( i=0; i<m; i++ ) { t[i] = b[i] - scale*(v[*diag++ + shift])*x[i]; } 990 PLogFlops(2*m); 991 992 993 /* t = (E + L)^{-1}t */ 994 ts = t + shift; /* shifted by one for index start of a or a->j*/ 995 diag = a->diag; 996 for ( i=0; i<m; i++ ) { 997 d = fshift + a->a[diag[i]+shift]; 998 n = diag[i] - a->i[i]; 999 PLogFlops(2*n-1); 1000 idx = a->j + a->i[i] + shift; 1001 v = a->a + a->i[i] + shift; 1002 sum = t[i]; 1003 SPARSEDENSEMDOT(sum,ts,v,idx,n); 1004 t[i] = omega*(sum/d); 1005 } 1006 1007 /* x = x + t */ 1008 for ( i=0; i<m; i++ ) { x[i] += t[i]; } 1009 PetscFree(t); 1010 ierr = VecRestoreArray(xx,&x); CHKERRQ(ierr); 1011 ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr); 1012 PetscFunctionReturn(0); 1013 } 1014 if (flag & SOR_ZERO_INITIAL_GUESS) { 1015 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1016 for ( i=0; i<m; i++ ) { 1017 d = fshift + a->a[diag[i]+shift]; 1018 n = diag[i] - a->i[i]; 1019 PLogFlops(2*n-1); 1020 idx = a->j + a->i[i] + shift; 1021 v = a->a + a->i[i] + shift; 1022 sum = b[i]; 1023 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1024 x[i] = omega*(sum/d); 1025 } 1026 xb = x; 1027 } else xb = b; 1028 if ((flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) && 1029 (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP)) { 1030 for ( i=0; i<m; i++ ) { 1031 x[i] *= a->a[diag[i]+shift]; 1032 } 1033 PLogFlops(m); 1034 } 1035 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1036 for ( i=m-1; i>=0; i-- ) { 1037 d = fshift + a->a[diag[i] + shift]; 1038 n = a->i[i+1] - diag[i] - 1; 1039 PLogFlops(2*n-1); 1040 idx = a->j + diag[i] + (!shift); 1041 v = a->a + diag[i] + (!shift); 1042 sum = xb[i]; 1043 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1044 x[i] = omega*(sum/d); 1045 } 1046 } 1047 its--; 1048 } 1049 while (its--) { 1050 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1051 for ( i=0; i<m; i++ ) { 1052 d = fshift + a->a[diag[i]+shift]; 1053 n = a->i[i+1] - a->i[i]; 1054 PLogFlops(2*n-1); 1055 idx = a->j + a->i[i] + shift; 1056 v = a->a + a->i[i] + shift; 1057 sum = b[i]; 1058 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1059 x[i] = (1. - omega)*x[i] + omega*(sum + a->a[diag[i]+shift]*x[i])/d; 1060 } 1061 } 1062 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1063 for ( i=m-1; i>=0; i-- ) { 1064 d = fshift + a->a[diag[i] + shift]; 1065 n = a->i[i+1] - a->i[i]; 1066 PLogFlops(2*n-1); 1067 idx = a->j + a->i[i] + shift; 1068 v = a->a + a->i[i] + shift; 1069 sum = b[i]; 1070 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1071 x[i] = (1. - omega)*x[i] + omega*(sum + a->a[diag[i]+shift]*x[i])/d; 1072 } 1073 } 1074 } 1075 ierr = VecRestoreArray(xx,&x); CHKERRQ(ierr); 1076 ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr); 1077 PetscFunctionReturn(0); 1078 } 1079 1080 #undef __FUNC__ 1081 #define __FUNC__ "MatGetInfo_SeqAIJ" 1082 int MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1083 { 1084 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1085 1086 PetscFunctionBegin; 1087 info->rows_global = (double)a->m; 1088 info->columns_global = (double)a->n; 1089 info->rows_local = (double)a->m; 1090 info->columns_local = (double)a->n; 1091 info->block_size = 1.0; 1092 info->nz_allocated = (double)a->maxnz; 1093 info->nz_used = (double)a->nz; 1094 info->nz_unneeded = (double)(a->maxnz - a->nz); 1095 info->assemblies = (double)A->num_ass; 1096 info->mallocs = (double)a->reallocs; 1097 info->memory = A->mem; 1098 if (A->factor) { 1099 info->fill_ratio_given = A->info.fill_ratio_given; 1100 info->fill_ratio_needed = A->info.fill_ratio_needed; 1101 info->factor_mallocs = A->info.factor_mallocs; 1102 } else { 1103 info->fill_ratio_given = 0; 1104 info->fill_ratio_needed = 0; 1105 info->factor_mallocs = 0; 1106 } 1107 PetscFunctionReturn(0); 1108 } 1109 1110 extern int MatLUFactorSymbolic_SeqAIJ(Mat,IS,IS,double,Mat*); 1111 extern int MatLUFactorNumeric_SeqAIJ(Mat,Mat*); 1112 extern int MatLUFactor_SeqAIJ(Mat,IS,IS,double); 1113 extern int MatSolve_SeqAIJ(Mat,Vec,Vec); 1114 extern int MatSolveAdd_SeqAIJ(Mat,Vec,Vec,Vec); 1115 extern int MatSolveTrans_SeqAIJ(Mat,Vec,Vec); 1116 extern int MatSolveTransAdd_SeqAIJ(Mat,Vec,Vec,Vec); 1117 1118 #undef __FUNC__ 1119 #define __FUNC__ "MatZeroRows_SeqAIJ" 1120 int MatZeroRows_SeqAIJ(Mat A,IS is,Scalar *diag) 1121 { 1122 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1123 int i,ierr,N, *rows,m = a->m - 1,shift = a->indexshift; 1124 1125 PetscFunctionBegin; 1126 ierr = ISGetSize(is,&N); CHKERRQ(ierr); 1127 ierr = ISGetIndices(is,&rows); CHKERRQ(ierr); 1128 if (diag) { 1129 for ( i=0; i<N; i++ ) { 1130 if (rows[i] < 0 || rows[i] > m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"row out of range"); 1131 if (a->ilen[rows[i]] > 0) { /* in case row was completely empty */ 1132 a->ilen[rows[i]] = 1; 1133 a->a[a->i[rows[i]]+shift] = *diag; 1134 a->j[a->i[rows[i]]+shift] = rows[i]+shift; 1135 } else { 1136 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],diag,INSERT_VALUES);CHKERRQ(ierr); 1137 } 1138 } 1139 } else { 1140 for ( i=0; i<N; i++ ) { 1141 if (rows[i] < 0 || rows[i] > m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"row out of range"); 1142 a->ilen[rows[i]] = 0; 1143 } 1144 } 1145 ISRestoreIndices(is,&rows); 1146 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1147 PetscFunctionReturn(0); 1148 } 1149 1150 #undef __FUNC__ 1151 #define __FUNC__ "MatGetSize_SeqAIJ" 1152 int MatGetSize_SeqAIJ(Mat A,int *m,int *n) 1153 { 1154 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1155 1156 PetscFunctionBegin; 1157 if (m) *m = a->m; 1158 if (n) *n = a->n; 1159 PetscFunctionReturn(0); 1160 } 1161 1162 #undef __FUNC__ 1163 #define __FUNC__ "MatGetOwnershipRange_SeqAIJ" 1164 int MatGetOwnershipRange_SeqAIJ(Mat A,int *m,int *n) 1165 { 1166 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1167 1168 PetscFunctionBegin; 1169 *m = 0; *n = a->m; 1170 PetscFunctionReturn(0); 1171 } 1172 1173 #undef __FUNC__ 1174 #define __FUNC__ "MatGetRow_SeqAIJ" 1175 int MatGetRow_SeqAIJ(Mat A,int row,int *nz,int **idx,Scalar **v) 1176 { 1177 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1178 int *itmp,i,shift = a->indexshift; 1179 1180 PetscFunctionBegin; 1181 if (row < 0 || row >= a->m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Row out of range"); 1182 1183 *nz = a->i[row+1] - a->i[row]; 1184 if (v) *v = a->a + a->i[row] + shift; 1185 if (idx) { 1186 itmp = a->j + a->i[row] + shift; 1187 if (*nz && shift) { 1188 *idx = (int *) PetscMalloc( (*nz)*sizeof(int) ); CHKPTRQ(*idx); 1189 for ( i=0; i<(*nz); i++ ) {(*idx)[i] = itmp[i] + shift;} 1190 } else if (*nz) { 1191 *idx = itmp; 1192 } 1193 else *idx = 0; 1194 } 1195 PetscFunctionReturn(0); 1196 } 1197 1198 #undef __FUNC__ 1199 #define __FUNC__ "MatRestoreRow_SeqAIJ" 1200 int MatRestoreRow_SeqAIJ(Mat A,int row,int *nz,int **idx,Scalar **v) 1201 { 1202 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1203 1204 PetscFunctionBegin; 1205 if (idx) {if (*idx && a->indexshift) PetscFree(*idx);} 1206 PetscFunctionReturn(0); 1207 } 1208 1209 #undef __FUNC__ 1210 #define __FUNC__ "MatNorm_SeqAIJ" 1211 int MatNorm_SeqAIJ(Mat A,NormType type,double *norm) 1212 { 1213 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1214 Scalar *v = a->a; 1215 double sum = 0.0; 1216 int i, j,shift = a->indexshift; 1217 1218 PetscFunctionBegin; 1219 if (type == NORM_FROBENIUS) { 1220 for (i=0; i<a->nz; i++ ) { 1221 #if defined(USE_PETSC_COMPLEX) 1222 sum += PetscReal(PetscConj(*v)*(*v)); v++; 1223 #else 1224 sum += (*v)*(*v); v++; 1225 #endif 1226 } 1227 *norm = sqrt(sum); 1228 } else if (type == NORM_1) { 1229 double *tmp; 1230 int *jj = a->j; 1231 tmp = (double *) PetscMalloc( (a->n+1)*sizeof(double) ); CHKPTRQ(tmp); 1232 PetscMemzero(tmp,a->n*sizeof(double)); 1233 *norm = 0.0; 1234 for ( j=0; j<a->nz; j++ ) { 1235 tmp[*jj++ + shift] += PetscAbsScalar(*v); v++; 1236 } 1237 for ( j=0; j<a->n; j++ ) { 1238 if (tmp[j] > *norm) *norm = tmp[j]; 1239 } 1240 PetscFree(tmp); 1241 } else if (type == NORM_INFINITY) { 1242 *norm = 0.0; 1243 for ( j=0; j<a->m; j++ ) { 1244 v = a->a + a->i[j] + shift; 1245 sum = 0.0; 1246 for ( i=0; i<a->i[j+1]-a->i[j]; i++ ) { 1247 sum += PetscAbsScalar(*v); v++; 1248 } 1249 if (sum > *norm) *norm = sum; 1250 } 1251 } else { 1252 SETERRQ(PETSC_ERR_SUP,0,"No support for two norm"); 1253 } 1254 PetscFunctionReturn(0); 1255 } 1256 1257 #undef __FUNC__ 1258 #define __FUNC__ "MatTranspose_SeqAIJ" 1259 int MatTranspose_SeqAIJ(Mat A,Mat *B) 1260 { 1261 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1262 Mat C; 1263 int i, ierr, *aj = a->j, *ai = a->i, m = a->m, len, *col; 1264 int shift = a->indexshift; 1265 Scalar *array = a->a; 1266 1267 PetscFunctionBegin; 1268 if (B == PETSC_NULL && m != a->n) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Square matrix only for in-place"); 1269 col = (int *) PetscMalloc((1+a->n)*sizeof(int)); CHKPTRQ(col); 1270 PetscMemzero(col,(1+a->n)*sizeof(int)); 1271 if (shift) { 1272 for ( i=0; i<ai[m]-1; i++ ) aj[i] -= 1; 1273 } 1274 for ( i=0; i<ai[m]+shift; i++ ) col[aj[i]] += 1; 1275 ierr = MatCreateSeqAIJ(A->comm,a->n,m,0,col,&C); CHKERRQ(ierr); 1276 PetscFree(col); 1277 for ( i=0; i<m; i++ ) { 1278 len = ai[i+1]-ai[i]; 1279 ierr = MatSetValues(C,len,aj,1,&i,array,INSERT_VALUES); CHKERRQ(ierr); 1280 array += len; aj += len; 1281 } 1282 if (shift) { 1283 for ( i=0; i<ai[m]-1; i++ ) aj[i] += 1; 1284 } 1285 1286 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1287 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1288 1289 if (B != PETSC_NULL) { 1290 *B = C; 1291 } else { 1292 PetscOps *Abops; 1293 MatOps Aops; 1294 1295 /* This isn't really an in-place transpose */ 1296 PetscFree(a->a); 1297 if (!a->singlemalloc) {PetscFree(a->i); PetscFree(a->j);} 1298 if (a->diag) PetscFree(a->diag); 1299 if (a->ilen) PetscFree(a->ilen); 1300 if (a->imax) PetscFree(a->imax); 1301 if (a->solve_work) PetscFree(a->solve_work); 1302 if (a->inode.size) PetscFree(a->inode.size); 1303 PetscFree(a); 1304 1305 1306 ierr = MapDestroy(A->rmap);CHKERRQ(ierr); 1307 ierr = MapDestroy(A->cmap);CHKERRQ(ierr); 1308 1309 /* 1310 This is horrible, horrible code. We need to keep the 1311 the bops and ops Structures, copy everything from C 1312 including the function pointers.. 1313 */ 1314 PetscMemcpy(A->ops,C->ops,sizeof(struct _MatOps)); 1315 PetscMemcpy(A->bops,C->bops,sizeof(PetscOps)); 1316 Abops = A->bops; 1317 Aops = A->ops; 1318 PetscMemcpy(A,C,sizeof(struct _p_Mat)); 1319 A->bops = Abops; 1320 A->ops = Aops; 1321 A->qlist = 0; 1322 1323 PetscHeaderDestroy(C); 1324 } 1325 PetscFunctionReturn(0); 1326 } 1327 1328 #undef __FUNC__ 1329 #define __FUNC__ "MatDiagonalScale_SeqAIJ" 1330 int MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 1331 { 1332 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1333 Scalar *l,*r,x,*v; 1334 int ierr,i,j,m = a->m, n = a->n, M, nz = a->nz, *jj,shift = a->indexshift; 1335 1336 PetscFunctionBegin; 1337 if (ll) { 1338 /* The local size is used so that VecMPI can be passed to this routine 1339 by MatDiagonalScale_MPIAIJ */ 1340 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 1341 if (m != a->m) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Left scaling vector wrong length"); 1342 ierr = VecGetArray(ll,&l); CHKERRQ(ierr); 1343 v = a->a; 1344 for ( i=0; i<m; i++ ) { 1345 x = l[i]; 1346 M = a->i[i+1] - a->i[i]; 1347 for ( j=0; j<M; j++ ) { (*v++) *= x;} 1348 } 1349 ierr = VecRestoreArray(ll,&l); CHKERRQ(ierr); 1350 PLogFlops(nz); 1351 } 1352 if (rr) { 1353 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 1354 if (n != a->n) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Right scaling vector wrong length"); 1355 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1356 v = a->a; jj = a->j; 1357 for ( i=0; i<nz; i++ ) { 1358 (*v++) *= r[*jj++ + shift]; 1359 } 1360 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1361 PLogFlops(nz); 1362 } 1363 PetscFunctionReturn(0); 1364 } 1365 1366 #undef __FUNC__ 1367 #define __FUNC__ "MatGetSubMatrix_SeqAIJ" 1368 int MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,int csize,MatGetSubMatrixCall scall,Mat *B) 1369 { 1370 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data,*c; 1371 int *smap, i, k, kstart, kend, ierr, oldcols = a->n,*lens; 1372 int row,mat_i,*mat_j,tcol,first,step,*mat_ilen; 1373 register int sum,lensi; 1374 int *irow, *icol, nrows, ncols, shift = a->indexshift,*ssmap; 1375 int *starts,*j_new,*i_new,*aj = a->j, *ai = a->i,ii,*ailen = a->ilen; 1376 Scalar *a_new,*mat_a; 1377 Mat C; 1378 PetscTruth stride; 1379 1380 PetscFunctionBegin; 1381 ierr = ISSorted(isrow,(PetscTruth*)&i);CHKERRQ(ierr); 1382 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"ISrow is not sorted"); 1383 ierr = ISSorted(iscol,(PetscTruth*)&i);CHKERRQ(ierr); 1384 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"IScol is not sorted"); 1385 1386 ierr = ISGetIndices(isrow,&irow); CHKERRQ(ierr); 1387 ierr = ISGetSize(isrow,&nrows); CHKERRQ(ierr); 1388 ierr = ISGetSize(iscol,&ncols); CHKERRQ(ierr); 1389 1390 ierr = ISStrideGetInfo(iscol,&first,&step); CHKERRQ(ierr); 1391 ierr = ISStride(iscol,&stride); CHKERRQ(ierr); 1392 if (stride && step == 1) { 1393 /* special case of contiguous rows */ 1394 lens = (int *) PetscMalloc((ncols+nrows+1)*sizeof(int)); CHKPTRQ(lens); 1395 starts = lens + ncols; 1396 /* loop over new rows determining lens and starting points */ 1397 for (i=0; i<nrows; i++) { 1398 kstart = ai[irow[i]]+shift; 1399 kend = kstart + ailen[irow[i]]; 1400 for ( k=kstart; k<kend; k++ ) { 1401 if (aj[k]+shift >= first) { 1402 starts[i] = k; 1403 break; 1404 } 1405 } 1406 sum = 0; 1407 while (k < kend) { 1408 if (aj[k++]+shift >= first+ncols) break; 1409 sum++; 1410 } 1411 lens[i] = sum; 1412 } 1413 /* create submatrix */ 1414 if (scall == MAT_REUSE_MATRIX) { 1415 int n_cols,n_rows; 1416 ierr = MatGetSize(*B,&n_rows,&n_cols); CHKERRQ(ierr); 1417 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Reused submatrix wrong size"); 1418 ierr = MatZeroEntries(*B); CHKERRQ(ierr); 1419 C = *B; 1420 } else { 1421 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1422 } 1423 c = (Mat_SeqAIJ*) C->data; 1424 1425 /* loop over rows inserting into submatrix */ 1426 a_new = c->a; 1427 j_new = c->j; 1428 i_new = c->i; 1429 i_new[0] = -shift; 1430 for (i=0; i<nrows; i++) { 1431 ii = starts[i]; 1432 lensi = lens[i]; 1433 for ( k=0; k<lensi; k++ ) { 1434 *j_new++ = aj[ii+k] - first; 1435 } 1436 PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(Scalar)); 1437 a_new += lensi; 1438 i_new[i+1] = i_new[i] + lensi; 1439 c->ilen[i] = lensi; 1440 } 1441 PetscFree(lens); 1442 } else { 1443 ierr = ISGetIndices(iscol,&icol); CHKERRQ(ierr); 1444 smap = (int *) PetscMalloc((1+oldcols)*sizeof(int)); CHKPTRQ(smap); 1445 ssmap = smap + shift; 1446 lens = (int *) PetscMalloc((1+nrows)*sizeof(int)); CHKPTRQ(lens); 1447 PetscMemzero(smap,oldcols*sizeof(int)); 1448 for ( i=0; i<ncols; i++ ) smap[icol[i]] = i+1; 1449 /* determine lens of each row */ 1450 for (i=0; i<nrows; i++) { 1451 kstart = ai[irow[i]]+shift; 1452 kend = kstart + a->ilen[irow[i]]; 1453 lens[i] = 0; 1454 for ( k=kstart; k<kend; k++ ) { 1455 if (ssmap[aj[k]]) { 1456 lens[i]++; 1457 } 1458 } 1459 } 1460 /* Create and fill new matrix */ 1461 if (scall == MAT_REUSE_MATRIX) { 1462 c = (Mat_SeqAIJ *)((*B)->data); 1463 if (c->m != nrows || c->n != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Cannot reuse matrix. wrong size"); 1464 if (PetscMemcmp(c->ilen,lens, c->m *sizeof(int))) { 1465 SETERRQ(PETSC_ERR_ARG_SIZ,0,"Cannot reuse matrix. wrong no of nonzeros"); 1466 } 1467 PetscMemzero(c->ilen,c->m*sizeof(int)); 1468 C = *B; 1469 } else { 1470 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1471 } 1472 c = (Mat_SeqAIJ *)(C->data); 1473 for (i=0; i<nrows; i++) { 1474 row = irow[i]; 1475 kstart = ai[row]+shift; 1476 kend = kstart + a->ilen[row]; 1477 mat_i = c->i[i]+shift; 1478 mat_j = c->j + mat_i; 1479 mat_a = c->a + mat_i; 1480 mat_ilen = c->ilen + i; 1481 for ( k=kstart; k<kend; k++ ) { 1482 if ((tcol=ssmap[a->j[k]])) { 1483 *mat_j++ = tcol - (!shift); 1484 *mat_a++ = a->a[k]; 1485 (*mat_ilen)++; 1486 1487 } 1488 } 1489 } 1490 /* Free work space */ 1491 ierr = ISRestoreIndices(iscol,&icol); CHKERRQ(ierr); 1492 PetscFree(smap); PetscFree(lens); 1493 } 1494 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1495 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1496 1497 ierr = ISRestoreIndices(isrow,&irow); CHKERRQ(ierr); 1498 *B = C; 1499 PetscFunctionReturn(0); 1500 } 1501 1502 /* 1503 note: This can only work for identity for row and col. It would 1504 be good to check this and otherwise generate an error. 1505 */ 1506 #undef __FUNC__ 1507 #define __FUNC__ "MatILUFactor_SeqAIJ" 1508 int MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,double efill,int fill) 1509 { 1510 Mat_SeqAIJ *a = (Mat_SeqAIJ *) inA->data; 1511 int ierr; 1512 Mat outA; 1513 1514 PetscFunctionBegin; 1515 if (fill != 0) SETERRQ(PETSC_ERR_SUP,0,"Only fill=0 supported"); 1516 1517 outA = inA; 1518 inA->factor = FACTOR_LU; 1519 a->row = row; 1520 a->col = col; 1521 1522 /* Create the invert permutation so that it can be used in MatLUFactorNumeric() */ 1523 ierr = ISInvertPermutation(col,&(a->icol)); CHKERRQ(ierr); 1524 PLogObjectParent(inA,a->icol); 1525 1526 if (!a->solve_work) { /* this matrix may have been factored before */ 1527 a->solve_work = (Scalar *) PetscMalloc( (a->m+1)*sizeof(Scalar));CHKPTRQ(a->solve_work); 1528 } 1529 1530 if (!a->diag) { 1531 ierr = MatMarkDiag_SeqAIJ(inA); CHKERRQ(ierr); 1532 } 1533 ierr = MatLUFactorNumeric_SeqAIJ(inA,&outA); CHKERRQ(ierr); 1534 PetscFunctionReturn(0); 1535 } 1536 1537 #include "pinclude/blaslapack.h" 1538 #undef __FUNC__ 1539 #define __FUNC__ "MatScale_SeqAIJ" 1540 int MatScale_SeqAIJ(Scalar *alpha,Mat inA) 1541 { 1542 Mat_SeqAIJ *a = (Mat_SeqAIJ *) inA->data; 1543 int one = 1; 1544 1545 PetscFunctionBegin; 1546 BLscal_( &a->nz, alpha, a->a, &one ); 1547 PLogFlops(a->nz); 1548 PetscFunctionReturn(0); 1549 } 1550 1551 #undef __FUNC__ 1552 #define __FUNC__ "MatGetSubMatrices_SeqAIJ" 1553 int MatGetSubMatrices_SeqAIJ(Mat A,int n, IS *irow,IS *icol,MatGetSubMatrixCall scall,Mat **B) 1554 { 1555 int ierr,i; 1556 1557 PetscFunctionBegin; 1558 if (scall == MAT_INITIAL_MATRIX) { 1559 *B = (Mat *) PetscMalloc( (n+1)*sizeof(Mat) ); CHKPTRQ(*B); 1560 } 1561 1562 for ( i=0; i<n; i++ ) { 1563 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1564 } 1565 PetscFunctionReturn(0); 1566 } 1567 1568 #undef __FUNC__ 1569 #define __FUNC__ "MatGetBlockSize_SeqAIJ" 1570 int MatGetBlockSize_SeqAIJ(Mat A, int *bs) 1571 { 1572 PetscFunctionBegin; 1573 *bs = 1; 1574 PetscFunctionReturn(0); 1575 } 1576 1577 #undef __FUNC__ 1578 #define __FUNC__ "MatIncreaseOverlap_SeqAIJ" 1579 int MatIncreaseOverlap_SeqAIJ(Mat A, int is_max, IS *is, int ov) 1580 { 1581 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1582 int shift, row, i,j,k,l,m,n, *idx,ierr, *nidx, isz, val; 1583 int start, end, *ai, *aj; 1584 BT table; 1585 1586 PetscFunctionBegin; 1587 shift = a->indexshift; 1588 m = a->m; 1589 ai = a->i; 1590 aj = a->j+shift; 1591 1592 if (ov < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"illegal overlap value used"); 1593 1594 nidx = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(nidx); 1595 ierr = BTCreate(m,table); CHKERRQ(ierr); 1596 1597 for ( i=0; i<is_max; i++ ) { 1598 /* Initialize the two local arrays */ 1599 isz = 0; 1600 BTMemzero(m,table); 1601 1602 /* Extract the indices, assume there can be duplicate entries */ 1603 ierr = ISGetIndices(is[i],&idx); CHKERRQ(ierr); 1604 ierr = ISGetSize(is[i],&n); CHKERRQ(ierr); 1605 1606 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 1607 for ( j=0; j<n ; ++j){ 1608 if(!BTLookupSet(table, idx[j])) { nidx[isz++] = idx[j];} 1609 } 1610 ierr = ISRestoreIndices(is[i],&idx); CHKERRQ(ierr); 1611 ierr = ISDestroy(is[i]); CHKERRQ(ierr); 1612 1613 k = 0; 1614 for ( j=0; j<ov; j++){ /* for each overlap */ 1615 n = isz; 1616 for ( ; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 1617 row = nidx[k]; 1618 start = ai[row]; 1619 end = ai[row+1]; 1620 for ( l = start; l<end ; l++){ 1621 val = aj[l] + shift; 1622 if (!BTLookupSet(table,val)) {nidx[isz++] = val;} 1623 } 1624 } 1625 } 1626 ierr = ISCreateGeneral(PETSC_COMM_SELF, isz, nidx, (is+i)); CHKERRQ(ierr); 1627 } 1628 BTDestroy(table); 1629 PetscFree(nidx); 1630 PetscFunctionReturn(0); 1631 } 1632 1633 /* -------------------------------------------------------------- */ 1634 #undef __FUNC__ 1635 #define __FUNC__ "MatPermute_SeqAIJ" 1636 int MatPermute_SeqAIJ(Mat A, IS rowp, IS colp, Mat *B) 1637 { 1638 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1639 Scalar *vwork; 1640 int i, ierr, nz, m = a->m, n = a->n, *cwork; 1641 int *row,*col,*cnew,j,*lens; 1642 IS icolp,irowp; 1643 1644 PetscFunctionBegin; 1645 ierr = ISInvertPermutation(rowp,&irowp); CHKERRQ(ierr); 1646 ierr = ISGetIndices(irowp,&row); CHKERRQ(ierr); 1647 ierr = ISInvertPermutation(colp,&icolp); CHKERRQ(ierr); 1648 ierr = ISGetIndices(icolp,&col); CHKERRQ(ierr); 1649 1650 /* determine lengths of permuted rows */ 1651 lens = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(lens); 1652 for (i=0; i<m; i++ ) { 1653 lens[row[i]] = a->i[i+1] - a->i[i]; 1654 } 1655 ierr = MatCreateSeqAIJ(A->comm,m,n,0,lens,B);CHKERRQ(ierr); 1656 PetscFree(lens); 1657 1658 cnew = (int *) PetscMalloc( n*sizeof(int) ); CHKPTRQ(cnew); 1659 for (i=0; i<m; i++) { 1660 ierr = MatGetRow(A,i,&nz,&cwork,&vwork); CHKERRQ(ierr); 1661 for (j=0; j<nz; j++ ) { cnew[j] = col[cwork[j]];} 1662 ierr = MatSetValues(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES); CHKERRQ(ierr); 1663 ierr = MatRestoreRow(A,i,&nz,&cwork,&vwork); CHKERRQ(ierr); 1664 } 1665 PetscFree(cnew); 1666 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1667 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 1668 ierr = ISRestoreIndices(irowp,&row); CHKERRQ(ierr); 1669 ierr = ISRestoreIndices(icolp,&col); CHKERRQ(ierr); 1670 ierr = ISDestroy(irowp); CHKERRQ(ierr); 1671 ierr = ISDestroy(icolp); CHKERRQ(ierr); 1672 PetscFunctionReturn(0); 1673 } 1674 1675 #undef __FUNC__ 1676 #define __FUNC__ "MatPrintHelp_SeqAIJ" 1677 int MatPrintHelp_SeqAIJ(Mat A) 1678 { 1679 static int called = 0; 1680 MPI_Comm comm = A->comm; 1681 1682 PetscFunctionBegin; 1683 if (called) {PetscFunctionReturn(0);} else called = 1; 1684 (*PetscHelpPrintf)(comm," Options for MATSEQAIJ and MATMPIAIJ matrix formats (the defaults):\n"); 1685 (*PetscHelpPrintf)(comm," -mat_lu_pivotthreshold <threshold>: Set pivoting threshold\n"); 1686 (*PetscHelpPrintf)(comm," -mat_aij_oneindex: internal indices begin at 1 instead of the default 0.\n"); 1687 (*PetscHelpPrintf)(comm," -mat_aij_no_inode: Do not use inodes\n"); 1688 (*PetscHelpPrintf)(comm," -mat_aij_inode_limit <limit>: Set inode limit (max limit=5)\n"); 1689 #if defined(HAVE_ESSL) 1690 (*PetscHelpPrintf)(comm," -mat_aij_essl: Use IBM sparse LU factorization and solve.\n"); 1691 #endif 1692 PetscFunctionReturn(0); 1693 } 1694 extern int MatEqual_SeqAIJ(Mat A,Mat B, PetscTruth* flg); 1695 extern int MatFDColoringCreate_SeqAIJ(Mat,ISColoring,MatFDColoring); 1696 extern int MatColoringPatch_SeqAIJ(Mat,int,int *,ISColoring *); 1697 1698 #undef __FUNC__ 1699 #define __FUNC__ "MatCopy_SeqAIJ" 1700 int MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 1701 { 1702 int ierr; 1703 1704 PetscFunctionBegin; 1705 if (str == SAME_NONZERO_PATTERN && B->type == MATSEQAIJ) { 1706 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 1707 Mat_SeqAIJ *b = (Mat_SeqAIJ *) B->data; 1708 1709 if (a->nonew != 1) SETERRQ(1,1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 1710 if (b->nonew != 1) SETERRQ(1,1,"Must call MatSetOption(B,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 1711 if (a->i[a->m]+a->indexshift != b->i[b->m]+a->indexshift) { 1712 SETERRQ(1,1,"Number of nonzeros in two matrices are different"); 1713 } 1714 PetscMemcpy(b->a,a->a,(a->i[a->m]+a->indexshift)*sizeof(Scalar)); 1715 } else { 1716 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1717 } 1718 PetscFunctionReturn(0); 1719 } 1720 1721 /* -------------------------------------------------------------------*/ 1722 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 1723 MatGetRow_SeqAIJ, 1724 MatRestoreRow_SeqAIJ, 1725 MatMult_SeqAIJ, 1726 MatMultAdd_SeqAIJ, 1727 MatMultTrans_SeqAIJ, 1728 MatMultTransAdd_SeqAIJ, 1729 MatSolve_SeqAIJ, 1730 MatSolveAdd_SeqAIJ, 1731 MatSolveTrans_SeqAIJ, 1732 MatSolveTransAdd_SeqAIJ, 1733 MatLUFactor_SeqAIJ, 1734 0, 1735 MatRelax_SeqAIJ, 1736 MatTranspose_SeqAIJ, 1737 MatGetInfo_SeqAIJ, 1738 MatEqual_SeqAIJ, 1739 MatGetDiagonal_SeqAIJ, 1740 MatDiagonalScale_SeqAIJ, 1741 MatNorm_SeqAIJ, 1742 0, 1743 MatAssemblyEnd_SeqAIJ, 1744 MatCompress_SeqAIJ, 1745 MatSetOption_SeqAIJ, 1746 MatZeroEntries_SeqAIJ, 1747 MatZeroRows_SeqAIJ, 1748 MatLUFactorSymbolic_SeqAIJ, 1749 MatLUFactorNumeric_SeqAIJ, 1750 0, 1751 0, 1752 MatGetSize_SeqAIJ, 1753 MatGetSize_SeqAIJ, 1754 MatGetOwnershipRange_SeqAIJ, 1755 MatILUFactorSymbolic_SeqAIJ, 1756 0, 1757 0, 1758 0, 1759 MatDuplicate_SeqAIJ, 1760 0, 1761 0, 1762 MatILUFactor_SeqAIJ, 1763 0, 1764 0, 1765 MatGetSubMatrices_SeqAIJ, 1766 MatIncreaseOverlap_SeqAIJ, 1767 MatGetValues_SeqAIJ, 1768 MatCopy_SeqAIJ, 1769 MatPrintHelp_SeqAIJ, 1770 MatScale_SeqAIJ, 1771 0, 1772 0, 1773 MatILUDTFactor_SeqAIJ, 1774 MatGetBlockSize_SeqAIJ, 1775 MatGetRowIJ_SeqAIJ, 1776 MatRestoreRowIJ_SeqAIJ, 1777 MatGetColumnIJ_SeqAIJ, 1778 MatRestoreColumnIJ_SeqAIJ, 1779 MatFDColoringCreate_SeqAIJ, 1780 MatColoringPatch_SeqAIJ, 1781 0, 1782 MatPermute_SeqAIJ, 1783 0, 1784 0, 1785 0, 1786 0, 1787 MatGetMaps_Petsc}; 1788 1789 extern int MatUseSuperLU_SeqAIJ(Mat); 1790 extern int MatUseEssl_SeqAIJ(Mat); 1791 extern int MatUseDXML_SeqAIJ(Mat); 1792 1793 #undef __FUNC__ 1794 #define __FUNC__ "MatSeqAIJSetColumnIndices_SeqAIJ" 1795 int MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,int *indices) 1796 { 1797 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 1798 int i,nz,n; 1799 1800 PetscFunctionBegin; 1801 if (aij->indexshift) SETERRQ(1,1,"No support with 1 based indexing"); 1802 1803 nz = aij->maxnz; 1804 n = aij->n; 1805 for (i=0; i<nz; i++) { 1806 aij->j[i] = indices[i]; 1807 } 1808 aij->nz = nz; 1809 for ( i=0; i<n; i++ ) { 1810 aij->ilen[i] = aij->imax[i]; 1811 } 1812 1813 PetscFunctionReturn(0); 1814 } 1815 1816 #undef __FUNC__ 1817 #define __FUNC__ "MatSeqAIJSetColumnIndices" 1818 /*@ 1819 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 1820 in the matrix. 1821 1822 Input Parameters: 1823 + mat - the SeqAIJ matrix 1824 - indices - the column indices 1825 1826 Notes: 1827 This can be called if you have precomputed the nonzero structure of the 1828 matrix and want to provide it to the matrix object to improve the performance 1829 of the MatSetValues() operation. 1830 1831 You MUST have set the correct numbers of nonzeros per row in the call to 1832 MatCreateSeqAIJ(). 1833 1834 MUST be called before any calls to MatSetValues(); 1835 1836 @*/ 1837 int MatSeqAIJSetColumnIndices(Mat mat,int *indices) 1838 { 1839 int ierr,(*f)(Mat,int *); 1840 1841 PetscFunctionBegin; 1842 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1843 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqAIJSetColumnIndices_C",(void **)&f); 1844 CHKERRQ(ierr); 1845 if (f) { 1846 ierr = (*f)(mat,indices);CHKERRQ(ierr); 1847 } else { 1848 SETERRQ(1,1,"Wrong type of matrix to set column indices"); 1849 } 1850 PetscFunctionReturn(0); 1851 } 1852 1853 /* ----------------------------------------------------------------------------------------*/ 1854 1855 #undef __FUNC__ 1856 #define __FUNC__ "MatStoreValues_SeqAIJ" 1857 int MatStoreValues_SeqAIJ(Mat mat) 1858 { 1859 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 1860 int nz = aij->i[aij->m]+aij->indexshift; 1861 1862 PetscFunctionBegin; 1863 if (aij->nonew != 1) { 1864 SETERRQ(1,1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 1865 } 1866 1867 /* allocate space for values if not already there */ 1868 if (!aij->saved_values) { 1869 aij->saved_values = (Scalar *) PetscMalloc(nz*sizeof(Scalar));CHKPTRQ(aij->saved_values); 1870 } 1871 1872 /* copy values over */ 1873 PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(Scalar)); 1874 PetscFunctionReturn(0); 1875 } 1876 1877 #undef __FUNC__ 1878 #define __FUNC__ "MatStoreValues" 1879 /*@ 1880 MatStoreValues - Stashes a copy of the matrix values; this allows, for 1881 example, reuse of the linear part of a Jacobian, while recomputing the 1882 nonlinear portion. 1883 1884 Collect on Mat 1885 1886 Input Parameters: 1887 . mat - the matrix (currently on AIJ matrices support this option) 1888 1889 Common Usage, with SNESSolve(): 1890 $ Create Jacobian matrix 1891 $ Set linear terms into matrix 1892 $ Apply boundary conditions to matrix, at this time matrix must have 1893 $ final nonzero structure (i.e. setting the nonlinear terms and applying 1894 $ boundary conditions again will not change the nonzero structure 1895 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 1896 $ ierr = MatStoreValues(mat); 1897 $ Call SNESSetJacobian() with matrix 1898 $ In your Jacobian routine 1899 $ ierr = MatRetrieveValues(mat); 1900 $ Set nonlinear terms in matrix 1901 1902 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 1903 $ // build linear portion of Jacobian 1904 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 1905 $ ierr = MatStoreValues(mat); 1906 $ loop over nonlinear iterations 1907 $ ierr = MatRetrieveValues(mat); 1908 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 1909 $ // call MatAssemblyBegin/End() on matrix 1910 $ Solve linear system with Jacobian 1911 $ endloop 1912 1913 Notes: 1914 Matrix must already be assemblied before calling this routine 1915 Must set the matrix option MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); before 1916 calling this routine. 1917 1918 .seealso: MatRetrieveValues() 1919 1920 @*/ 1921 int MatStoreValues(Mat mat) 1922 { 1923 int ierr,(*f)(Mat); 1924 1925 PetscFunctionBegin; 1926 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1927 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1928 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 1929 1930 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatStoreValues_C",(void **)&f); 1931 CHKERRQ(ierr); 1932 if (f) { 1933 ierr = (*f)(mat);CHKERRQ(ierr); 1934 } else { 1935 SETERRQ(1,1,"Wrong type of matrix to store values"); 1936 } 1937 PetscFunctionReturn(0); 1938 } 1939 1940 #undef __FUNC__ 1941 #define __FUNC__ "MatRetrieveValues_SeqAIJ" 1942 int MatRetrieveValues_SeqAIJ(Mat mat) 1943 { 1944 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 1945 int nz = aij->i[aij->m]+aij->indexshift; 1946 1947 PetscFunctionBegin; 1948 if (aij->nonew != 1) { 1949 SETERRQ(1,1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 1950 } 1951 if (!aij->saved_values) { 1952 SETERRQ(1,1,"Must call MatStoreValues(A);first"); 1953 } 1954 1955 /* copy values over */ 1956 PetscMemcpy(aij->a, aij->saved_values,nz*sizeof(Scalar)); 1957 PetscFunctionReturn(0); 1958 } 1959 1960 #undef __FUNC__ 1961 #define __FUNC__ "MatRetrieveValues" 1962 /*@ 1963 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 1964 example, reuse of the linear part of a Jacobian, while recomputing the 1965 nonlinear portion. 1966 1967 Collect on Mat 1968 1969 Input Parameters: 1970 . mat - the matrix (currently on AIJ matrices support this option) 1971 1972 .seealso: MatStoreValues() 1973 1974 @*/ 1975 int MatRetrieveValues(Mat mat) 1976 { 1977 int ierr,(*f)(Mat); 1978 1979 PetscFunctionBegin; 1980 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1981 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 1982 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 1983 1984 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatRetrieveValues_C",(void **)&f); 1985 CHKERRQ(ierr); 1986 if (f) { 1987 ierr = (*f)(mat);CHKERRQ(ierr); 1988 } else { 1989 SETERRQ(1,1,"Wrong type of matrix to retrieve values"); 1990 } 1991 PetscFunctionReturn(0); 1992 } 1993 1994 /* --------------------------------------------------------------------------------*/ 1995 1996 #undef __FUNC__ 1997 #define __FUNC__ "MatCreateSeqAIJ" 1998 /*@C 1999 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 2000 (the default parallel PETSc format). For good matrix assembly performance 2001 the user should preallocate the matrix storage by setting the parameter nz 2002 (or the array nzz). By setting these parameters accurately, performance 2003 during matrix assembly can be increased by more than a factor of 50. 2004 2005 Collective on MPI_Comm 2006 2007 Input Parameters: 2008 + comm - MPI communicator, set to PETSC_COMM_SELF 2009 . m - number of rows 2010 . n - number of columns 2011 . nz - number of nonzeros per row (same for all rows) 2012 - nzz - array containing the number of nonzeros in the various rows 2013 (possibly different for each row) or PETSC_NULL 2014 2015 Output Parameter: 2016 . A - the matrix 2017 2018 Notes: 2019 The AIJ format (also called the Yale sparse matrix format or 2020 compressed row storage), is fully compatible with standard Fortran 77 2021 storage. That is, the stored row and column indices can begin at 2022 either one (as in Fortran) or zero. See the users' manual for details. 2023 2024 Specify the preallocated storage with either nz or nnz (not both). 2025 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2026 allocation. For large problems you MUST preallocate memory or you 2027 will get TERRIBLE performance, see the users' manual chapter on matrices. 2028 2029 By default, this format uses inodes (identical nodes) when possible, to 2030 improve numerical efficiency of matrix-vector products and solves. We 2031 search for consecutive rows with the same nonzero structure, thereby 2032 reusing matrix information to achieve increased efficiency. 2033 2034 Options Database Keys: 2035 + -mat_aij_no_inode - Do not use inodes 2036 . -mat_aij_inode_limit <limit> - Sets inode limit (max limit=5) 2037 - -mat_aij_oneindex - Internally use indexing starting at 1 2038 rather than 0. Note that when calling MatSetValues(), 2039 the user still MUST index entries starting at 0! 2040 2041 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices() 2042 @*/ 2043 int MatCreateSeqAIJ(MPI_Comm comm,int m,int n,int nz,int *nnz, Mat *A) 2044 { 2045 Mat B; 2046 Mat_SeqAIJ *b; 2047 int i, len, ierr, flg,size; 2048 2049 PetscFunctionBegin; 2050 MPI_Comm_size(comm,&size); 2051 if (size > 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Comm must be of size 1"); 2052 2053 *A = 0; 2054 PetscHeaderCreate(B,_p_Mat,struct _MatOps,MAT_COOKIE,MATSEQAIJ,comm,MatDestroy,MatView); 2055 PLogObjectCreate(B); 2056 B->data = (void *) (b = PetscNew(Mat_SeqAIJ)); CHKPTRQ(b); 2057 PetscMemzero(b,sizeof(Mat_SeqAIJ)); 2058 PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps)); 2059 B->ops->destroy = MatDestroy_SeqAIJ; 2060 B->ops->view = MatView_SeqAIJ; 2061 B->factor = 0; 2062 B->lupivotthreshold = 1.0; 2063 B->mapping = 0; 2064 ierr = OptionsGetDouble(PETSC_NULL,"-mat_lu_pivotthreshold",&B->lupivotthreshold,&flg);CHKERRQ(ierr); 2065 b->ilu_preserve_row_sums = PETSC_FALSE; 2066 ierr = OptionsHasName(PETSC_NULL,"-pc_ilu_preserve_row_sums",(int*)&b->ilu_preserve_row_sums);CHKERRQ(ierr); 2067 b->row = 0; 2068 b->col = 0; 2069 b->icol = 0; 2070 b->indexshift = 0; 2071 b->reallocs = 0; 2072 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_oneindex", &flg); CHKERRQ(ierr); 2073 if (flg) b->indexshift = -1; 2074 2075 b->m = m; B->m = m; B->M = m; 2076 b->n = n; B->n = n; B->N = n; 2077 2078 ierr = MapCreateMPI(comm,m,m,&B->rmap);CHKERRQ(ierr); 2079 ierr = MapCreateMPI(comm,n,n,&B->cmap);CHKERRQ(ierr); 2080 2081 b->imax = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(b->imax); 2082 if (nnz == PETSC_NULL) { 2083 if (nz == PETSC_DEFAULT) nz = 10; 2084 else if (nz <= 0) nz = 1; 2085 for ( i=0; i<m; i++ ) b->imax[i] = nz; 2086 nz = nz*m; 2087 } else { 2088 nz = 0; 2089 for ( i=0; i<m; i++ ) {b->imax[i] = nnz[i]; nz += nnz[i];} 2090 } 2091 2092 /* allocate the matrix space */ 2093 len = nz*(sizeof(int) + sizeof(Scalar)) + (b->m+1)*sizeof(int); 2094 b->a = (Scalar *) PetscMalloc( len ); CHKPTRQ(b->a); 2095 b->j = (int *) (b->a + nz); 2096 PetscMemzero(b->j,nz*sizeof(int)); 2097 b->i = b->j + nz; 2098 b->singlemalloc = 1; 2099 2100 b->i[0] = -b->indexshift; 2101 for (i=1; i<m+1; i++) { 2102 b->i[i] = b->i[i-1] + b->imax[i-1]; 2103 } 2104 2105 /* b->ilen will count nonzeros in each row so far. */ 2106 b->ilen = (int *) PetscMalloc((m+1)*sizeof(int)); 2107 PLogObjectMemory(B,len+2*(m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2108 for ( i=0; i<b->m; i++ ) { b->ilen[i] = 0;} 2109 2110 b->nz = 0; 2111 b->maxnz = nz; 2112 b->sorted = 0; 2113 b->roworiented = 1; 2114 b->nonew = 0; 2115 b->diag = 0; 2116 b->solve_work = 0; 2117 b->spptr = 0; 2118 b->inode.node_count = 0; 2119 b->inode.size = 0; 2120 b->inode.limit = 5; 2121 b->inode.max_limit = 5; 2122 b->saved_values = 0; 2123 B->info.nz_unneeded = (double)b->maxnz; 2124 2125 *A = B; 2126 2127 /* SuperLU is not currently supported through PETSc */ 2128 #if defined(HAVE_SUPERLU) 2129 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_superlu", &flg); CHKERRQ(ierr); 2130 if (flg) { ierr = MatUseSuperLU_SeqAIJ(B); CHKERRQ(ierr); } 2131 #endif 2132 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_essl", &flg); CHKERRQ(ierr); 2133 if (flg) { ierr = MatUseEssl_SeqAIJ(B); CHKERRQ(ierr); } 2134 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_dxml", &flg); CHKERRQ(ierr); 2135 if (flg) { 2136 if (!b->indexshift) SETERRQ( PETSC_ERR_LIB,0,"need -mat_aij_oneindex with -mat_aij_dxml"); 2137 ierr = MatUseDXML_SeqAIJ(B); CHKERRQ(ierr); 2138 } 2139 ierr = OptionsHasName(PETSC_NULL,"-help", &flg); CHKERRQ(ierr); 2140 if (flg) {ierr = MatPrintHelp(B); CHKERRQ(ierr); } 2141 2142 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 2143 "MatSeqAIJSetColumnIndices_SeqAIJ", 2144 (void*)MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 2145 ierr = PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C", 2146 "MatStoreValues_SeqAIJ", 2147 (void*)MatStoreValues_SeqAIJ);CHKERRQ(ierr); 2148 ierr = PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C", 2149 "MatRetrieveValues_SeqAIJ", 2150 (void*)MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 2151 PetscFunctionReturn(0); 2152 } 2153 2154 #undef __FUNC__ 2155 #define __FUNC__ "MatDuplicate_SeqAIJ" 2156 int MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 2157 { 2158 Mat C; 2159 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ *) A->data; 2160 int i,len, m = a->m,shift = a->indexshift,ierr; 2161 2162 PetscFunctionBegin; 2163 *B = 0; 2164 PetscHeaderCreate(C,_p_Mat,struct _MatOps,MAT_COOKIE,MATSEQAIJ,A->comm,MatDestroy,MatView); 2165 PLogObjectCreate(C); 2166 C->data = (void *) (c = PetscNew(Mat_SeqAIJ)); CHKPTRQ(c); 2167 PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps)); 2168 C->ops->destroy = MatDestroy_SeqAIJ; 2169 C->ops->view = MatView_SeqAIJ; 2170 C->factor = A->factor; 2171 c->row = 0; 2172 c->col = 0; 2173 c->icol = 0; 2174 c->indexshift = shift; 2175 C->assembled = PETSC_TRUE; 2176 2177 c->m = C->m = a->m; 2178 c->n = C->n = a->n; 2179 C->M = a->m; 2180 C->N = a->n; 2181 2182 c->imax = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(c->imax); 2183 c->ilen = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(c->ilen); 2184 for ( i=0; i<m; i++ ) { 2185 c->imax[i] = a->imax[i]; 2186 c->ilen[i] = a->ilen[i]; 2187 } 2188 2189 /* allocate the matrix space */ 2190 c->singlemalloc = 1; 2191 len = (m+1)*sizeof(int)+(a->i[m])*(sizeof(Scalar)+sizeof(int)); 2192 c->a = (Scalar *) PetscMalloc( len ); CHKPTRQ(c->a); 2193 c->j = (int *) (c->a + a->i[m] + shift); 2194 c->i = c->j + a->i[m] + shift; 2195 PetscMemcpy(c->i,a->i,(m+1)*sizeof(int)); 2196 if (m > 0) { 2197 PetscMemcpy(c->j,a->j,(a->i[m]+shift)*sizeof(int)); 2198 if (cpvalues == MAT_COPY_VALUES) { 2199 PetscMemcpy(c->a,a->a,(a->i[m]+shift)*sizeof(Scalar)); 2200 } else { 2201 PetscMemzero(c->a,(a->i[m]+shift)*sizeof(Scalar)); 2202 } 2203 } 2204 2205 PLogObjectMemory(C,len+2*(m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2206 c->sorted = a->sorted; 2207 c->roworiented = a->roworiented; 2208 c->nonew = a->nonew; 2209 c->ilu_preserve_row_sums = a->ilu_preserve_row_sums; 2210 c->saved_values = 0; 2211 2212 if (a->diag) { 2213 c->diag = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(c->diag); 2214 PLogObjectMemory(C,(m+1)*sizeof(int)); 2215 for ( i=0; i<m; i++ ) { 2216 c->diag[i] = a->diag[i]; 2217 } 2218 } else c->diag = 0; 2219 c->inode.limit = a->inode.limit; 2220 c->inode.max_limit = a->inode.max_limit; 2221 if (a->inode.size){ 2222 c->inode.size = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(c->inode.size); 2223 c->inode.node_count = a->inode.node_count; 2224 PetscMemcpy( c->inode.size, a->inode.size, (m+1)*sizeof(int)); 2225 } else { 2226 c->inode.size = 0; 2227 c->inode.node_count = 0; 2228 } 2229 c->nz = a->nz; 2230 c->maxnz = a->maxnz; 2231 c->solve_work = 0; 2232 c->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 2233 2234 *B = C; 2235 ierr = PetscObjectComposeFunction((PetscObject)C,"MatSeqAIJSetColumnIndices_C", 2236 "MatSeqAIJSetColumnIndices_SeqAIJ", 2237 (void*)MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 2238 PetscFunctionReturn(0); 2239 } 2240 2241 #undef __FUNC__ 2242 #define __FUNC__ "MatLoad_SeqAIJ" 2243 int MatLoad_SeqAIJ(Viewer viewer,MatType type,Mat *A) 2244 { 2245 Mat_SeqAIJ *a; 2246 Mat B; 2247 int i, nz, ierr, fd, header[4],size,*rowlengths = 0,M,N,shift; 2248 MPI_Comm comm; 2249 2250 PetscFunctionBegin; 2251 PetscObjectGetComm((PetscObject) viewer,&comm); 2252 MPI_Comm_size(comm,&size); 2253 if (size > 1) SETERRQ(PETSC_ERR_ARG_SIZ,0,"view must have one processor"); 2254 ierr = ViewerBinaryGetDescriptor(viewer,&fd); CHKERRQ(ierr); 2255 ierr = PetscBinaryRead(fd,header,4,PETSC_INT); CHKERRQ(ierr); 2256 if (header[0] != MAT_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,0,"not matrix object in file"); 2257 M = header[1]; N = header[2]; nz = header[3]; 2258 2259 if (nz < 0) { 2260 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,1,"Matrix stored in special format on disk, cannot load as SeqAIJ"); 2261 } 2262 2263 /* read in row lengths */ 2264 rowlengths = (int*) PetscMalloc( M*sizeof(int) ); CHKPTRQ(rowlengths); 2265 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT); CHKERRQ(ierr); 2266 2267 /* create our matrix */ 2268 ierr = MatCreateSeqAIJ(comm,M,N,0,rowlengths,A); CHKERRQ(ierr); 2269 B = *A; 2270 a = (Mat_SeqAIJ *) B->data; 2271 shift = a->indexshift; 2272 2273 /* read in column indices and adjust for Fortran indexing*/ 2274 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT); CHKERRQ(ierr); 2275 if (shift) { 2276 for ( i=0; i<nz; i++ ) { 2277 a->j[i] += 1; 2278 } 2279 } 2280 2281 /* read in nonzero values */ 2282 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR); CHKERRQ(ierr); 2283 2284 /* set matrix "i" values */ 2285 a->i[0] = -shift; 2286 for ( i=1; i<= M; i++ ) { 2287 a->i[i] = a->i[i-1] + rowlengths[i-1]; 2288 a->ilen[i-1] = rowlengths[i-1]; 2289 } 2290 PetscFree(rowlengths); 2291 2292 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 2293 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 2294 PetscFunctionReturn(0); 2295 } 2296 2297 #undef __FUNC__ 2298 #define __FUNC__ "MatEqual_SeqAIJ" 2299 int MatEqual_SeqAIJ(Mat A,Mat B, PetscTruth* flg) 2300 { 2301 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)B->data; 2302 2303 PetscFunctionBegin; 2304 if (B->type !=MATSEQAIJ)SETERRQ(PETSC_ERR_ARG_INCOMP,0,"Matrices must be same type"); 2305 2306 /* If the matrix dimensions are not equal, or no of nonzeros or shift */ 2307 if ((a->m != b->m ) || (a->n !=b->n) ||( a->nz != b->nz)|| 2308 (a->indexshift != b->indexshift)) { 2309 *flg = PETSC_FALSE; PetscFunctionReturn(0); 2310 } 2311 2312 /* if the a->i are the same */ 2313 if (PetscMemcmp(a->i,b->i,(a->m+1)*sizeof(int))) { 2314 *flg = PETSC_FALSE; PetscFunctionReturn(0); 2315 } 2316 2317 /* if a->j are the same */ 2318 if (PetscMemcmp(a->j, b->j, (a->nz)*sizeof(int))) { 2319 *flg = PETSC_FALSE; PetscFunctionReturn(0); 2320 } 2321 2322 /* if a->a are the same */ 2323 if (PetscMemcmp(a->a, b->a, (a->nz)*sizeof(Scalar))) { 2324 *flg = PETSC_FALSE; PetscFunctionReturn(0); 2325 } 2326 *flg = PETSC_TRUE; 2327 PetscFunctionReturn(0); 2328 2329 } 2330