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