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