1 #ifdef PETSC_RCS_HEADER 2 static char vcid[] = "$Id: aij.c,v 1.302 1999/03/08 22:38:38 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 (a->ssor) PetscFree(a->ssor); 706 #if defined(USE_PETSC_LOG) 707 PLogObjectState((PetscObject)A,"Rows=%d, Cols=%d, NZ=%d",a->m,a->n,a->nz); 708 #endif 709 PetscFree(a->a); 710 if (!a->singlemalloc) { PetscFree(a->i); PetscFree(a->j);} 711 if (a->diag) PetscFree(a->diag); 712 if (a->ilen) PetscFree(a->ilen); 713 if (a->imax) PetscFree(a->imax); 714 if (a->solve_work) PetscFree(a->solve_work); 715 if (a->inode.size) PetscFree(a->inode.size); 716 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} 717 if (a->saved_values) PetscFree(a->saved_values); 718 PetscFree(a); 719 720 PLogObjectDestroy(A); 721 PetscHeaderDestroy(A); 722 PetscFunctionReturn(0); 723 } 724 725 #undef __FUNC__ 726 #define __FUNC__ "MatCompress_SeqAIJ" 727 int MatCompress_SeqAIJ(Mat A) 728 { 729 PetscFunctionBegin; 730 PetscFunctionReturn(0); 731 } 732 733 #undef __FUNC__ 734 #define __FUNC__ "MatSetOption_SeqAIJ" 735 int MatSetOption_SeqAIJ(Mat A,MatOption op) 736 { 737 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 738 739 PetscFunctionBegin; 740 if (op == MAT_ROW_ORIENTED) a->roworiented = 1; 741 else if (op == MAT_COLUMN_ORIENTED) a->roworiented = 0; 742 else if (op == MAT_COLUMNS_SORTED) a->sorted = 1; 743 else if (op == MAT_COLUMNS_UNSORTED) a->sorted = 0; 744 else if (op == MAT_NO_NEW_NONZERO_LOCATIONS) a->nonew = 1; 745 else if (op == MAT_NEW_NONZERO_LOCATION_ERR) a->nonew = -1; 746 else if (op == MAT_NEW_NONZERO_ALLOCATION_ERR) a->nonew = -2; 747 else if (op == MAT_YES_NEW_NONZERO_LOCATIONS) a->nonew = 0; 748 else if (op == MAT_ROWS_SORTED || 749 op == MAT_ROWS_UNSORTED || 750 op == MAT_SYMMETRIC || 751 op == MAT_STRUCTURALLY_SYMMETRIC || 752 op == MAT_YES_NEW_DIAGONALS || 753 op == MAT_IGNORE_OFF_PROC_ENTRIES|| 754 op == MAT_USE_HASH_TABLE) 755 PLogInfo(A,"MatSetOption_SeqAIJ:Option ignored\n"); 756 else if (op == MAT_NO_NEW_DIAGONALS) { 757 SETERRQ(PETSC_ERR_SUP,0,"MAT_NO_NEW_DIAGONALS"); 758 } else if (op == MAT_INODE_LIMIT_1) a->inode.limit = 1; 759 else if (op == MAT_INODE_LIMIT_2) a->inode.limit = 2; 760 else if (op == MAT_INODE_LIMIT_3) a->inode.limit = 3; 761 else if (op == MAT_INODE_LIMIT_4) a->inode.limit = 4; 762 else if (op == MAT_INODE_LIMIT_5) a->inode.limit = 5; 763 else SETERRQ(PETSC_ERR_SUP,0,"unknown option"); 764 PetscFunctionReturn(0); 765 } 766 767 #undef __FUNC__ 768 #define __FUNC__ "MatGetDiagonal_SeqAIJ" 769 int MatGetDiagonal_SeqAIJ(Mat A,Vec v) 770 { 771 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 772 int i,j, n,shift = a->indexshift,ierr; 773 Scalar *x, zero = 0.0; 774 775 PetscFunctionBegin; 776 ierr = VecSet(&zero,v);CHKERRQ(ierr); 777 ierr = VecGetArray(v,&x);;CHKERRQ(ierr); 778 ierr = VecGetLocalSize(v,&n);;CHKERRQ(ierr); 779 if (n != a->m) SETERRQ(PETSC_ERR_ARG_SIZ,0,"Nonconforming matrix and vector"); 780 for ( i=0; i<a->m; i++ ) { 781 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 782 if (a->j[j]+shift == i) { 783 x[i] = a->a[j]; 784 break; 785 } 786 } 787 } 788 ierr = VecRestoreArray(v,&x);;CHKERRQ(ierr); 789 PetscFunctionReturn(0); 790 } 791 792 /* -------------------------------------------------------*/ 793 /* Should check that shapes of vectors and matrices match */ 794 /* -------------------------------------------------------*/ 795 #undef __FUNC__ 796 #define __FUNC__ "MatMultTrans_SeqAIJ" 797 int MatMultTrans_SeqAIJ(Mat A,Vec xx,Vec yy) 798 { 799 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 800 Scalar *x, *y, *v, alpha; 801 int ierr,m = a->m, n, i, *idx, shift = a->indexshift; 802 803 PetscFunctionBegin; 804 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 805 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 806 PetscMemzero(y,a->n*sizeof(Scalar)); 807 y = y + shift; /* shift for Fortran start by 1 indexing */ 808 for ( i=0; i<m; i++ ) { 809 idx = a->j + a->i[i] + shift; 810 v = a->a + a->i[i] + shift; 811 n = a->i[i+1] - a->i[i]; 812 alpha = x[i]; 813 while (n-->0) {y[*idx++] += alpha * *v++;} 814 } 815 PLogFlops(2*a->nz - a->n); 816 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 817 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 818 PetscFunctionReturn(0); 819 } 820 821 #undef __FUNC__ 822 #define __FUNC__ "MatMultTransAdd_SeqAIJ" 823 int MatMultTransAdd_SeqAIJ(Mat A,Vec xx,Vec zz,Vec yy) 824 { 825 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 826 Scalar *x, *y, *v, alpha; 827 int ierr,m = a->m, n, i, *idx,shift = a->indexshift; 828 829 PetscFunctionBegin; 830 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 831 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 832 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 833 y = y + shift; /* shift for Fortran start by 1 indexing */ 834 for ( i=0; i<m; i++ ) { 835 idx = a->j + a->i[i] + shift; 836 v = a->a + a->i[i] + shift; 837 n = a->i[i+1] - a->i[i]; 838 alpha = x[i]; 839 while (n-->0) {y[*idx++] += alpha * *v++;} 840 } 841 PLogFlops(2*a->nz); 842 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 843 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 844 PetscFunctionReturn(0); 845 } 846 847 #undef __FUNC__ 848 #define __FUNC__ "MatMult_SeqAIJ" 849 int MatMult_SeqAIJ(Mat A,Vec xx,Vec yy) 850 { 851 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 852 Scalar *x, *y, *v, sum; 853 int ierr,m = a->m, *idx, shift = a->indexshift,*ii; 854 #if !defined(USE_FORTRAN_KERNEL_MULTAIJ) 855 int n, i, jrow,j; 856 #endif 857 858 #if defined(HAVE_PRAGMA_DISJOINT) 859 #pragma disjoint(*x,*y,*v) 860 #endif 861 862 PetscFunctionBegin; 863 ierr = VecGetArray(xx,&x); CHKERRQ(ierr); 864 ierr = VecGetArray(yy,&y); CHKERRQ(ierr); 865 x = x + shift; /* shift for Fortran start by 1 indexing */ 866 idx = a->j; 867 v = a->a; 868 ii = a->i; 869 #if defined(USE_FORTRAN_KERNEL_MULTAIJ) 870 fortranmultaij_(&m,x,ii,idx+shift,v+shift,y); 871 #else 872 v += shift; /* shift for Fortran start by 1 indexing */ 873 idx += shift; 874 for ( i=0; i<m; i++ ) { 875 jrow = ii[i]; 876 n = ii[i+1] - jrow; 877 sum = 0.0; 878 for ( j=0; j<n; j++) { 879 sum += v[jrow]*x[idx[jrow]]; jrow++; 880 } 881 y[i] = sum; 882 } 883 #endif 884 PLogFlops(2*a->nz - m); 885 ierr = VecRestoreArray(xx,&x); CHKERRQ(ierr); 886 ierr = VecRestoreArray(yy,&y); CHKERRQ(ierr); 887 PetscFunctionReturn(0); 888 } 889 890 #undef __FUNC__ 891 #define __FUNC__ "MatMultAdd_SeqAIJ" 892 int MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 893 { 894 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 895 Scalar *x, *y, *z, *v, sum; 896 int ierr,m = a->m, *idx, shift = a->indexshift,*ii; 897 #if !defined(USE_FORTRAN_KERNEL_MULTADDAIJ) 898 int n,i,jrow,j; 899 #endif 900 901 PetscFunctionBegin; 902 ierr = VecGetArray(xx,&x); CHKERRQ(ierr); 903 ierr = VecGetArray(yy,&y); CHKERRQ(ierr); 904 if (zz != yy) { 905 ierr = VecGetArray(zz,&z); CHKERRQ(ierr); 906 } else { 907 z = y; 908 } 909 x = x + shift; /* shift for Fortran start by 1 indexing */ 910 idx = a->j; 911 v = a->a; 912 ii = a->i; 913 #if defined(USE_FORTRAN_KERNEL_MULTADDAIJ) 914 fortranmultaddaij_(&m,x,ii,idx+shift,v+shift,y,z); 915 #else 916 v += shift; /* shift for Fortran start by 1 indexing */ 917 idx += shift; 918 for ( i=0; i<m; i++ ) { 919 jrow = ii[i]; 920 n = ii[i+1] - jrow; 921 sum = y[i]; 922 for ( j=0; j<n; j++) { 923 sum += v[jrow]*x[idx[jrow]]; jrow++; 924 } 925 z[i] = sum; 926 } 927 #endif 928 PLogFlops(2*a->nz); 929 ierr = VecRestoreArray(xx,&x); CHKERRQ(ierr); 930 ierr = VecRestoreArray(yy,&y); CHKERRQ(ierr); 931 if (zz != yy) { 932 ierr = VecRestoreArray(zz,&z); CHKERRQ(ierr); 933 } 934 PetscFunctionReturn(0); 935 } 936 937 /* 938 Adds diagonal pointers to sparse matrix structure. 939 */ 940 #undef __FUNC__ 941 #define __FUNC__ "MatMarkDiag_SeqAIJ" 942 int MatMarkDiag_SeqAIJ(Mat A) 943 { 944 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 945 int i,j, *diag, m = a->m,shift = a->indexshift; 946 947 PetscFunctionBegin; 948 diag = (int *) PetscMalloc( (m+1)*sizeof(int)); CHKPTRQ(diag); 949 PLogObjectMemory(A,(m+1)*sizeof(int)); 950 for ( i=0; i<a->m; i++ ) { 951 diag[i] = a->i[i+1]; 952 for ( j=a->i[i]+shift; j<a->i[i+1]+shift; j++ ) { 953 if (a->j[j]+shift == i) { 954 diag[i] = j - shift; 955 break; 956 } 957 } 958 } 959 a->diag = diag; 960 PetscFunctionReturn(0); 961 } 962 963 /* 964 Checks for missing diagonals 965 */ 966 #undef __FUNC__ 967 #define __FUNC__ "MatMissingDiag_SeqAIJ" 968 int MatMissingDiag_SeqAIJ(Mat A) 969 { 970 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 971 int *diag = a->diag, *jj = a->j,i,shift = a->indexshift; 972 973 PetscFunctionBegin; 974 for ( i=0; i<a->m; i++ ) { 975 if (jj[diag[i]+shift] != i-shift) { 976 SETERRQ1(1,1,"Matrix is missing diagonal number %d",i); 977 } 978 } 979 PetscFunctionReturn(0); 980 } 981 982 #undef __FUNC__ 983 #define __FUNC__ "MatRelax_SeqAIJ" 984 int MatRelax_SeqAIJ(Mat A,Vec bb,double omega,MatSORType flag,double fshift,int its,Vec xx) 985 { 986 Mat_SeqAIJ *a = (Mat_SeqAIJ *) A->data; 987 Scalar *x, *b, *bs, d, *xs, sum, *v = a->a,*t,scale,*ts, *xb; 988 int ierr, *idx, *diag,n = a->n, m = a->m, i, shift = a->indexshift; 989 990 PetscFunctionBegin; 991 ierr = VecGetArray(xx,&x); CHKERRQ(ierr); 992 if (xx != bb) { 993 ierr = VecGetArray(bb,&b);CHKERRQ(ierr); 994 } else { 995 b = x; 996 } 997 998 if (!a->diag) {ierr = MatMarkDiag_SeqAIJ(A);CHKERRQ(ierr);} 999 diag = a->diag; 1000 xs = x + shift; /* shifted by one for index start of a or a->j*/ 1001 if (flag == SOR_APPLY_UPPER) { 1002 /* apply ( U + D/omega) to the vector */ 1003 bs = b + shift; 1004 for ( i=0; i<m; i++ ) { 1005 d = fshift + a->a[diag[i] + shift]; 1006 n = a->i[i+1] - diag[i] - 1; 1007 PLogFlops(2*n-1); 1008 idx = a->j + diag[i] + (!shift); 1009 v = a->a + diag[i] + (!shift); 1010 sum = b[i]*d/omega; 1011 SPARSEDENSEDOT(sum,bs,v,idx,n); 1012 x[i] = sum; 1013 } 1014 ierr = VecRestoreArray(xx,&x); CHKERRQ(ierr); 1015 if (bb != xx) {ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);} 1016 PetscFunctionReturn(0); 1017 } 1018 if (flag == SOR_APPLY_LOWER) { 1019 SETERRQ(PETSC_ERR_SUP,0,"SOR_APPLY_LOWER is not done"); 1020 } else if (0 && (flag & SOR_EISENSTAT) && omega == 1.0 && shift == 0 && fshift == 0.0) { 1021 /* Let A = L + U + D; where L is lower trianglar, 1022 U is upper triangular, E is diagonal; This routine applies 1023 1024 (L + E)^{-1} A (U + E)^{-1} 1025 1026 to a vector efficiently using Eisenstat's trick. This is for 1027 the case of SSOR preconditioner, so E is D/omega where omega 1028 is the relaxation factor; but in this special case omega == 1 1029 */ 1030 t = (Scalar *) PetscMalloc( m*sizeof(Scalar) ); CHKPTRQ(t); 1031 1032 /* x = (E + U)^{-1} b */ 1033 for ( i=m-1; i>=0; i-- ) { 1034 d = a->a[diag[i]]; 1035 n = a->i[i+1] - diag[i] - 1; 1036 PLogFlops(2*n-1); 1037 idx = a->j + diag[i] + 1; 1038 v = a->a + diag[i] + 1; 1039 sum = b[i]; 1040 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1041 x[i] = sum/d; 1042 } 1043 1044 /* t = b - (2*E - D)x */ 1045 v = a->a; 1046 for ( i=0; i<m; i++ ) { t[i] = b[i] - (v[*diag++])*x[i]; } 1047 PLogFlops(2*m); 1048 1049 1050 /* t = (E + L)^{-1}t */ 1051 diag = a->diag; 1052 for ( i=0; i<m; i++ ) { 1053 d = a->a[diag[i]]; 1054 n = diag[i] - a->i[i]; 1055 PLogFlops(2*n-1); 1056 idx = a->j + a->i[i]; 1057 v = a->a + a->i[i]; 1058 sum = t[i]; 1059 SPARSEDENSEMDOT(sum,t,v,idx,n); 1060 t[i] = sum/d; 1061 } 1062 1063 /* x = x + t */ 1064 for ( i=0; i<m; i++ ) { x[i] += t[i]; } 1065 PetscFree(t); 1066 ierr = VecRestoreArray(xx,&x); CHKERRQ(ierr); 1067 if (bb != xx) {ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);} 1068 PetscFunctionReturn(0); 1069 } else if (flag & SOR_EISENSTAT) { 1070 /* Let A = L + U + D; where L is lower trianglar, 1071 U is upper triangular, E is diagonal; This routine applies 1072 1073 (L + E)^{-1} A (U + E)^{-1} 1074 1075 to a vector efficiently using Eisenstat's trick. This is for 1076 the case of SSOR preconditioner, so E is D/omega where omega 1077 is the relaxation factor. 1078 */ 1079 t = (Scalar *) PetscMalloc( m*sizeof(Scalar) ); CHKPTRQ(t); 1080 scale = (2.0/omega) - 1.0; 1081 1082 /* x = (E + U)^{-1} b */ 1083 for ( i=m-1; i>=0; i-- ) { 1084 d = fshift + a->a[diag[i] + shift]; 1085 n = a->i[i+1] - diag[i] - 1; 1086 PLogFlops(2*n-1); 1087 idx = a->j + diag[i] + (!shift); 1088 v = a->a + diag[i] + (!shift); 1089 sum = b[i]; 1090 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1091 x[i] = omega*(sum/d); 1092 } 1093 1094 /* t = b - (2*E - D)x */ 1095 v = a->a; 1096 for ( i=0; i<m; i++ ) { t[i] = b[i] - scale*(v[*diag++ + shift])*x[i]; } 1097 PLogFlops(3*m); 1098 1099 1100 /* t = (E + L)^{-1}t */ 1101 ts = t + shift; /* shifted by one for index start of a or a->j*/ 1102 diag = a->diag; 1103 for ( i=0; i<m; i++ ) { 1104 d = fshift + a->a[diag[i]+shift]; 1105 n = diag[i] - a->i[i]; 1106 PLogFlops(2*n-1); 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 } 1113 1114 /* x = x + t */ 1115 for ( i=0; i<m; i++ ) { x[i] += t[i]; } 1116 PetscFree(t); 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 int MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,int *indices) 1902 { 1903 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 1904 int i,nz,n; 1905 1906 PetscFunctionBegin; 1907 if (aij->indexshift) SETERRQ(1,1,"No support with 1 based indexing"); 1908 1909 nz = aij->maxnz; 1910 n = aij->n; 1911 for (i=0; i<nz; i++) { 1912 aij->j[i] = indices[i]; 1913 } 1914 aij->nz = nz; 1915 for ( i=0; i<n; i++ ) { 1916 aij->ilen[i] = aij->imax[i]; 1917 } 1918 1919 PetscFunctionReturn(0); 1920 } 1921 EXTERN_C_END 1922 1923 #undef __FUNC__ 1924 #define __FUNC__ "MatSeqAIJSetColumnIndices" 1925 /*@ 1926 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 1927 in the matrix. 1928 1929 Input Parameters: 1930 + mat - the SeqAIJ matrix 1931 - indices - the column indices 1932 1933 Notes: 1934 This can be called if you have precomputed the nonzero structure of the 1935 matrix and want to provide it to the matrix object to improve the performance 1936 of the MatSetValues() operation. 1937 1938 You MUST have set the correct numbers of nonzeros per row in the call to 1939 MatCreateSeqAIJ(). 1940 1941 MUST be called before any calls to MatSetValues(); 1942 1943 @*/ 1944 int MatSeqAIJSetColumnIndices(Mat mat,int *indices) 1945 { 1946 int ierr,(*f)(Mat,int *); 1947 1948 PetscFunctionBegin; 1949 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1950 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqAIJSetColumnIndices_C",(void **)&f); 1951 CHKERRQ(ierr); 1952 if (f) { 1953 ierr = (*f)(mat,indices);CHKERRQ(ierr); 1954 } else { 1955 SETERRQ(1,1,"Wrong type of matrix to set column indices"); 1956 } 1957 PetscFunctionReturn(0); 1958 } 1959 1960 /* ----------------------------------------------------------------------------------------*/ 1961 1962 EXTERN_C_BEGIN 1963 #undef __FUNC__ 1964 #define __FUNC__ "MatStoreValues_SeqAIJ" 1965 int MatStoreValues_SeqAIJ(Mat mat) 1966 { 1967 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 1968 int nz = aij->i[aij->m]+aij->indexshift; 1969 1970 PetscFunctionBegin; 1971 if (aij->nonew != 1) { 1972 SETERRQ(1,1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 1973 } 1974 1975 /* allocate space for values if not already there */ 1976 if (!aij->saved_values) { 1977 aij->saved_values = (Scalar *) PetscMalloc(nz*sizeof(Scalar));CHKPTRQ(aij->saved_values); 1978 } 1979 1980 /* copy values over */ 1981 PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(Scalar)); 1982 PetscFunctionReturn(0); 1983 } 1984 EXTERN_C_END 1985 1986 #undef __FUNC__ 1987 #define __FUNC__ "MatStoreValues" 1988 /*@ 1989 MatStoreValues - Stashes a copy of the matrix values; this allows, for 1990 example, reuse of the linear part of a Jacobian, while recomputing the 1991 nonlinear portion. 1992 1993 Collect on Mat 1994 1995 Input Parameters: 1996 . mat - the matrix (currently on AIJ matrices support this option) 1997 1998 Common Usage, with SNESSolve(): 1999 $ Create Jacobian matrix 2000 $ Set linear terms into matrix 2001 $ Apply boundary conditions to matrix, at this time matrix must have 2002 $ final nonzero structure (i.e. setting the nonlinear terms and applying 2003 $ boundary conditions again will not change the nonzero structure 2004 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 2005 $ ierr = MatStoreValues(mat); 2006 $ Call SNESSetJacobian() with matrix 2007 $ In your Jacobian routine 2008 $ ierr = MatRetrieveValues(mat); 2009 $ Set nonlinear terms in matrix 2010 2011 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 2012 $ // build linear portion of Jacobian 2013 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 2014 $ ierr = MatStoreValues(mat); 2015 $ loop over nonlinear iterations 2016 $ ierr = MatRetrieveValues(mat); 2017 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 2018 $ // call MatAssemblyBegin/End() on matrix 2019 $ Solve linear system with Jacobian 2020 $ endloop 2021 2022 Notes: 2023 Matrix must already be assemblied before calling this routine 2024 Must set the matrix option MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); before 2025 calling this routine. 2026 2027 .seealso: MatRetrieveValues() 2028 2029 @*/ 2030 int MatStoreValues(Mat mat) 2031 { 2032 int ierr,(*f)(Mat); 2033 2034 PetscFunctionBegin; 2035 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2036 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 2037 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 2038 2039 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatStoreValues_C",(void **)&f); 2040 CHKERRQ(ierr); 2041 if (f) { 2042 ierr = (*f)(mat);CHKERRQ(ierr); 2043 } else { 2044 SETERRQ(1,1,"Wrong type of matrix to store values"); 2045 } 2046 PetscFunctionReturn(0); 2047 } 2048 2049 EXTERN_C_BEGIN 2050 #undef __FUNC__ 2051 #define __FUNC__ "MatRetrieveValues_SeqAIJ" 2052 int MatRetrieveValues_SeqAIJ(Mat mat) 2053 { 2054 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2055 int nz = aij->i[aij->m]+aij->indexshift; 2056 2057 PetscFunctionBegin; 2058 if (aij->nonew != 1) { 2059 SETERRQ(1,1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 2060 } 2061 if (!aij->saved_values) { 2062 SETERRQ(1,1,"Must call MatStoreValues(A);first"); 2063 } 2064 2065 /* copy values over */ 2066 PetscMemcpy(aij->a, aij->saved_values,nz*sizeof(Scalar)); 2067 PetscFunctionReturn(0); 2068 } 2069 EXTERN_C_END 2070 2071 #undef __FUNC__ 2072 #define __FUNC__ "MatRetrieveValues" 2073 /*@ 2074 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 2075 example, reuse of the linear part of a Jacobian, while recomputing the 2076 nonlinear portion. 2077 2078 Collect on Mat 2079 2080 Input Parameters: 2081 . mat - the matrix (currently on AIJ matrices support this option) 2082 2083 .seealso: MatStoreValues() 2084 2085 @*/ 2086 int MatRetrieveValues(Mat mat) 2087 { 2088 int ierr,(*f)(Mat); 2089 2090 PetscFunctionBegin; 2091 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2092 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for unassembled matrix"); 2093 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,0,"Not for factored matrix"); 2094 2095 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatRetrieveValues_C",(void **)&f); 2096 CHKERRQ(ierr); 2097 if (f) { 2098 ierr = (*f)(mat);CHKERRQ(ierr); 2099 } else { 2100 SETERRQ(1,1,"Wrong type of matrix to retrieve values"); 2101 } 2102 PetscFunctionReturn(0); 2103 } 2104 2105 /* --------------------------------------------------------------------------------*/ 2106 2107 #undef __FUNC__ 2108 #define __FUNC__ "MatCreateSeqAIJ" 2109 /*@C 2110 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 2111 (the default parallel PETSc format). For good matrix assembly performance 2112 the user should preallocate the matrix storage by setting the parameter nz 2113 (or the array nzz). By setting these parameters accurately, performance 2114 during matrix assembly can be increased by more than a factor of 50. 2115 2116 Collective on MPI_Comm 2117 2118 Input Parameters: 2119 + comm - MPI communicator, set to PETSC_COMM_SELF 2120 . m - number of rows 2121 . n - number of columns 2122 . nz - number of nonzeros per row (same for all rows) 2123 - nzz - array containing the number of nonzeros in the various rows 2124 (possibly different for each row) or PETSC_NULL 2125 2126 Output Parameter: 2127 . A - the matrix 2128 2129 Notes: 2130 The AIJ format (also called the Yale sparse matrix format or 2131 compressed row storage), is fully compatible with standard Fortran 77 2132 storage. That is, the stored row and column indices can begin at 2133 either one (as in Fortran) or zero. See the users' manual for details. 2134 2135 Specify the preallocated storage with either nz or nnz (not both). 2136 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2137 allocation. For large problems you MUST preallocate memory or you 2138 will get TERRIBLE performance, see the users' manual chapter on matrices. 2139 2140 By default, this format uses inodes (identical nodes) when possible, to 2141 improve numerical efficiency of matrix-vector products and solves. We 2142 search for consecutive rows with the same nonzero structure, thereby 2143 reusing matrix information to achieve increased efficiency. 2144 2145 Options Database Keys: 2146 + -mat_aij_no_inode - Do not use inodes 2147 . -mat_aij_inode_limit <limit> - Sets inode limit (max limit=5) 2148 - -mat_aij_oneindex - Internally use indexing starting at 1 2149 rather than 0. Note that when calling MatSetValues(), 2150 the user still MUST index entries starting at 0! 2151 2152 Level: intermediate 2153 2154 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices() 2155 @*/ 2156 int MatCreateSeqAIJ(MPI_Comm comm,int m,int n,int nz,int *nnz, Mat *A) 2157 { 2158 Mat B; 2159 Mat_SeqAIJ *b; 2160 int i, len, ierr, flg,size; 2161 2162 PetscFunctionBegin; 2163 MPI_Comm_size(comm,&size); 2164 if (size > 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,0,"Comm must be of size 1"); 2165 2166 *A = 0; 2167 PetscHeaderCreate(B,_p_Mat,struct _MatOps,MAT_COOKIE,MATSEQAIJ,"Mat",comm,MatDestroy,MatView); 2168 PLogObjectCreate(B); 2169 B->data = (void *) (b = PetscNew(Mat_SeqAIJ)); CHKPTRQ(b); 2170 PetscMemzero(b,sizeof(Mat_SeqAIJ)); 2171 PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps)); 2172 B->ops->destroy = MatDestroy_SeqAIJ; 2173 B->ops->view = MatView_SeqAIJ; 2174 B->factor = 0; 2175 B->lupivotthreshold = 1.0; 2176 B->mapping = 0; 2177 ierr = OptionsGetDouble(PETSC_NULL,"-mat_lu_pivotthreshold",&B->lupivotthreshold,&flg);CHKERRQ(ierr); 2178 b->ilu_preserve_row_sums = PETSC_FALSE; 2179 ierr = OptionsHasName(PETSC_NULL,"-pc_ilu_preserve_row_sums",(int*)&b->ilu_preserve_row_sums);CHKERRQ(ierr); 2180 b->row = 0; 2181 b->col = 0; 2182 b->icol = 0; 2183 b->indexshift = 0; 2184 b->reallocs = 0; 2185 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_oneindex", &flg); CHKERRQ(ierr); 2186 if (flg) b->indexshift = -1; 2187 2188 b->m = m; B->m = m; B->M = m; 2189 b->n = n; B->n = n; B->N = n; 2190 2191 ierr = MapCreateMPI(comm,m,m,&B->rmap);CHKERRQ(ierr); 2192 ierr = MapCreateMPI(comm,n,n,&B->cmap);CHKERRQ(ierr); 2193 2194 b->imax = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(b->imax); 2195 if (nnz == PETSC_NULL) { 2196 if (nz == PETSC_DEFAULT) nz = 10; 2197 else if (nz <= 0) nz = 1; 2198 for ( i=0; i<m; i++ ) b->imax[i] = nz; 2199 nz = nz*m; 2200 } else { 2201 nz = 0; 2202 for ( i=0; i<m; i++ ) {b->imax[i] = nnz[i]; nz += nnz[i];} 2203 } 2204 2205 /* allocate the matrix space */ 2206 len = nz*(sizeof(int) + sizeof(Scalar)) + (b->m+1)*sizeof(int); 2207 b->a = (Scalar *) PetscMalloc( len ); CHKPTRQ(b->a); 2208 b->j = (int *) (b->a + nz); 2209 PetscMemzero(b->j,nz*sizeof(int)); 2210 b->i = b->j + nz; 2211 b->singlemalloc = 1; 2212 2213 b->i[0] = -b->indexshift; 2214 for (i=1; i<m+1; i++) { 2215 b->i[i] = b->i[i-1] + b->imax[i-1]; 2216 } 2217 2218 /* b->ilen will count nonzeros in each row so far. */ 2219 b->ilen = (int *) PetscMalloc((m+1)*sizeof(int)); 2220 PLogObjectMemory(B,len+2*(m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2221 for ( i=0; i<b->m; i++ ) { b->ilen[i] = 0;} 2222 2223 b->nz = 0; 2224 b->maxnz = nz; 2225 b->sorted = 0; 2226 b->roworiented = 1; 2227 b->nonew = 0; 2228 b->diag = 0; 2229 b->solve_work = 0; 2230 b->spptr = 0; 2231 b->inode.node_count = 0; 2232 b->inode.size = 0; 2233 b->inode.limit = 5; 2234 b->inode.max_limit = 5; 2235 b->saved_values = 0; 2236 B->info.nz_unneeded = (double)b->maxnz; 2237 2238 *A = B; 2239 2240 /* SuperLU is not currently supported through PETSc */ 2241 #if defined(HAVE_SUPERLU) 2242 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_superlu", &flg); CHKERRQ(ierr); 2243 if (flg) { ierr = MatUseSuperLU_SeqAIJ(B); CHKERRQ(ierr); } 2244 #endif 2245 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_essl", &flg); CHKERRQ(ierr); 2246 if (flg) { ierr = MatUseEssl_SeqAIJ(B); CHKERRQ(ierr); } 2247 ierr = OptionsHasName(PETSC_NULL,"-mat_aij_dxml", &flg); CHKERRQ(ierr); 2248 if (flg) { 2249 if (!b->indexshift) SETERRQ( PETSC_ERR_LIB,0,"need -mat_aij_oneindex with -mat_aij_dxml"); 2250 ierr = MatUseDXML_SeqAIJ(B); CHKERRQ(ierr); 2251 } 2252 ierr = OptionsHasName(PETSC_NULL,"-help", &flg); CHKERRQ(ierr); 2253 if (flg) {ierr = MatPrintHelp(B); CHKERRQ(ierr); } 2254 2255 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 2256 "MatSeqAIJSetColumnIndices_SeqAIJ", 2257 (void*)MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 2258 ierr = PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C", 2259 "MatStoreValues_SeqAIJ", 2260 (void*)MatStoreValues_SeqAIJ);CHKERRQ(ierr); 2261 ierr = PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C", 2262 "MatRetrieveValues_SeqAIJ", 2263 (void*)MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 2264 PetscFunctionReturn(0); 2265 } 2266 2267 #undef __FUNC__ 2268 #define __FUNC__ "MatDuplicate_SeqAIJ" 2269 int MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 2270 { 2271 Mat C; 2272 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ *) A->data; 2273 int i,len, m = a->m,shift = a->indexshift,ierr; 2274 2275 PetscFunctionBegin; 2276 *B = 0; 2277 PetscHeaderCreate(C,_p_Mat,struct _MatOps,MAT_COOKIE,MATSEQAIJ,"Mat",A->comm,MatDestroy,MatView); 2278 PLogObjectCreate(C); 2279 C->data = (void *) (c = PetscNew(Mat_SeqAIJ)); CHKPTRQ(c); 2280 PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps)); 2281 C->ops->destroy = MatDestroy_SeqAIJ; 2282 C->ops->view = MatView_SeqAIJ; 2283 C->factor = A->factor; 2284 c->row = 0; 2285 c->col = 0; 2286 c->icol = 0; 2287 c->indexshift = shift; 2288 C->assembled = PETSC_TRUE; 2289 2290 c->m = C->m = a->m; 2291 c->n = C->n = a->n; 2292 C->M = a->m; 2293 C->N = a->n; 2294 2295 c->imax = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(c->imax); 2296 c->ilen = (int *) PetscMalloc((m+1)*sizeof(int)); CHKPTRQ(c->ilen); 2297 for ( i=0; i<m; i++ ) { 2298 c->imax[i] = a->imax[i]; 2299 c->ilen[i] = a->ilen[i]; 2300 } 2301 2302 /* allocate the matrix space */ 2303 c->singlemalloc = 1; 2304 len = (m+1)*sizeof(int)+(a->i[m])*(sizeof(Scalar)+sizeof(int)); 2305 c->a = (Scalar *) PetscMalloc( len ); CHKPTRQ(c->a); 2306 c->j = (int *) (c->a + a->i[m] + shift); 2307 c->i = c->j + a->i[m] + shift; 2308 PetscMemcpy(c->i,a->i,(m+1)*sizeof(int)); 2309 if (m > 0) { 2310 PetscMemcpy(c->j,a->j,(a->i[m]+shift)*sizeof(int)); 2311 if (cpvalues == MAT_COPY_VALUES) { 2312 PetscMemcpy(c->a,a->a,(a->i[m]+shift)*sizeof(Scalar)); 2313 } else { 2314 PetscMemzero(c->a,(a->i[m]+shift)*sizeof(Scalar)); 2315 } 2316 } 2317 2318 PLogObjectMemory(C,len+2*(m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2319 c->sorted = a->sorted; 2320 c->roworiented = a->roworiented; 2321 c->nonew = a->nonew; 2322 c->ilu_preserve_row_sums = a->ilu_preserve_row_sums; 2323 c->saved_values = 0; 2324 2325 if (a->diag) { 2326 c->diag = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(c->diag); 2327 PLogObjectMemory(C,(m+1)*sizeof(int)); 2328 for ( i=0; i<m; i++ ) { 2329 c->diag[i] = a->diag[i]; 2330 } 2331 } else c->diag = 0; 2332 c->inode.limit = a->inode.limit; 2333 c->inode.max_limit = a->inode.max_limit; 2334 if (a->inode.size){ 2335 c->inode.size = (int *) PetscMalloc( (m+1)*sizeof(int) ); CHKPTRQ(c->inode.size); 2336 c->inode.node_count = a->inode.node_count; 2337 PetscMemcpy( c->inode.size, a->inode.size, (m+1)*sizeof(int)); 2338 } else { 2339 c->inode.size = 0; 2340 c->inode.node_count = 0; 2341 } 2342 c->nz = a->nz; 2343 c->maxnz = a->maxnz; 2344 c->solve_work = 0; 2345 c->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 2346 2347 *B = C; 2348 ierr = FListDuplicate(A->qlist,&C->qlist);CHKERRQ(ierr); 2349 PetscFunctionReturn(0); 2350 } 2351 2352 #undef __FUNC__ 2353 #define __FUNC__ "MatLoad_SeqAIJ" 2354 int MatLoad_SeqAIJ(Viewer viewer,MatType type,Mat *A) 2355 { 2356 Mat_SeqAIJ *a; 2357 Mat B; 2358 int i, nz, ierr, fd, header[4],size,*rowlengths = 0,M,N,shift; 2359 MPI_Comm comm; 2360 2361 PetscFunctionBegin; 2362 ierr = PetscObjectGetComm((PetscObject) viewer,&comm);CHKERRQ(ierr); 2363 MPI_Comm_size(comm,&size); 2364 if (size > 1) SETERRQ(PETSC_ERR_ARG_SIZ,0,"view must have one processor"); 2365 ierr = ViewerBinaryGetDescriptor(viewer,&fd); CHKERRQ(ierr); 2366 ierr = PetscBinaryRead(fd,header,4,PETSC_INT); CHKERRQ(ierr); 2367 if (header[0] != MAT_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,0,"not matrix object in file"); 2368 M = header[1]; N = header[2]; nz = header[3]; 2369 2370 if (nz < 0) { 2371 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,1,"Matrix stored in special format on disk, cannot load as SeqAIJ"); 2372 } 2373 2374 /* read in row lengths */ 2375 rowlengths = (int*) PetscMalloc( M*sizeof(int) ); CHKPTRQ(rowlengths); 2376 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT); CHKERRQ(ierr); 2377 2378 /* create our matrix */ 2379 ierr = MatCreateSeqAIJ(comm,M,N,0,rowlengths,A); CHKERRQ(ierr); 2380 B = *A; 2381 a = (Mat_SeqAIJ *) B->data; 2382 shift = a->indexshift; 2383 2384 /* read in column indices and adjust for Fortran indexing*/ 2385 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT); CHKERRQ(ierr); 2386 if (shift) { 2387 for ( i=0; i<nz; i++ ) { 2388 a->j[i] += 1; 2389 } 2390 } 2391 2392 /* read in nonzero values */ 2393 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR); CHKERRQ(ierr); 2394 2395 /* set matrix "i" values */ 2396 a->i[0] = -shift; 2397 for ( i=1; i<= M; i++ ) { 2398 a->i[i] = a->i[i-1] + rowlengths[i-1]; 2399 a->ilen[i-1] = rowlengths[i-1]; 2400 } 2401 PetscFree(rowlengths); 2402 2403 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 2404 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr); 2405 PetscFunctionReturn(0); 2406 } 2407 2408 #undef __FUNC__ 2409 #define __FUNC__ "MatEqual_SeqAIJ" 2410 int MatEqual_SeqAIJ(Mat A,Mat B, PetscTruth* flg) 2411 { 2412 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data, *b = (Mat_SeqAIJ *)B->data; 2413 2414 PetscFunctionBegin; 2415 if (B->type !=MATSEQAIJ)SETERRQ(PETSC_ERR_ARG_INCOMP,0,"Matrices must be same type"); 2416 2417 /* If the matrix dimensions are not equal, or no of nonzeros or shift */ 2418 if ((a->m != b->m ) || (a->n !=b->n) ||( a->nz != b->nz)|| 2419 (a->indexshift != b->indexshift)) { 2420 *flg = PETSC_FALSE; PetscFunctionReturn(0); 2421 } 2422 2423 /* if the a->i are the same */ 2424 if (PetscMemcmp(a->i,b->i,(a->m+1)*sizeof(int))) { 2425 *flg = PETSC_FALSE; PetscFunctionReturn(0); 2426 } 2427 2428 /* if a->j are the same */ 2429 if (PetscMemcmp(a->j, b->j, (a->nz)*sizeof(int))) { 2430 *flg = PETSC_FALSE; PetscFunctionReturn(0); 2431 } 2432 2433 /* if a->a are the same */ 2434 if (PetscMemcmp(a->a, b->a, (a->nz)*sizeof(Scalar))) { 2435 *flg = PETSC_FALSE; PetscFunctionReturn(0); 2436 } 2437 *flg = PETSC_TRUE; 2438 PetscFunctionReturn(0); 2439 2440 } 2441