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