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