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