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