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