1 /*$Id: aij.c,v 1.365 2001/02/27 22:02:32 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__ "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 && !a->structurally_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 ierr = PetscMalloc((A->m+1)*sizeof(int),ia);CHKERRQ(ierr); 34 ierr = PetscMalloc((nz+1)*sizeof(int),ja);CHKERRQ(ierr); 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 ierr = PetscMalloc((A->m+1)*sizeof(int),ia);CHKERRQ(ierr); 41 ierr = PetscMalloc((nz+1)*sizeof(int),ja);CHKERRQ(ierr); 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__ "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 && !a->structurally_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__ "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 ierr = PetscMalloc((n+1)*sizeof(int),&collengths);CHKERRQ(ierr); 81 ierr = PetscMemzero(collengths,n*sizeof(int));CHKERRQ(ierr); 82 ierr = PetscMalloc((n+1)*sizeof(int),&cia);CHKERRQ(ierr); 83 ierr = PetscMalloc((nz+1)*sizeof(int),&cja);CHKERRQ(ierr); 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__ "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__ "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 ierr = PetscMalloc(len,&new_a);CHKERRQ(ierr); 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 PetscLogObjectMemory(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__ "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__ "MatView_SeqAIJ_Binary" 270 int MatView_SeqAIJ_Binary(Mat A,PetscViewer viewer) 271 { 272 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 273 int i,fd,*col_lens,ierr; 274 275 PetscFunctionBegin; 276 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 277 ierr = PetscMalloc((4+A->m)*sizeof(int),&col_lens);CHKERRQ(ierr); 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__ "MatView_SeqAIJ_ASCII" 306 int MatView_SeqAIJ_ASCII(Mat A,PetscViewer viewer) 307 { 308 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 309 int ierr,i,j,m = A->m,shift = a->indexshift; 310 char *name; 311 PetscViewerFormat format; 312 313 PetscFunctionBegin; 314 ierr = PetscObjectGetName((PetscObject)viewer,&name);CHKERRQ(ierr); 315 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 316 if (format == PETSC_VIEWER_ASCII_INFO_LONG || format == PETSC_VIEWER_ASCII_INFO) { 317 if (a->inode.size) { 318 ierr = PetscViewerASCIIPrintf(viewer,"using I-node routines: found %d nodes, limit used is %d\n",a->inode.node_count,a->inode.limit);CHKERRQ(ierr); 319 } else { 320 ierr = PetscViewerASCIIPrintf(viewer,"not using I-node routines\n");CHKERRQ(ierr); 321 } 322 } else if (format == PETSC_VIEWER_ASCII_MATLAB) { 323 int nofinalvalue = 0; 324 if ((a->i[m] == a->i[m-1]) || (a->j[a->nz-1] != A->n-!shift)) { 325 nofinalvalue = 1; 326 } 327 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 328 ierr = PetscViewerASCIIPrintf(viewer,"%% Size = %d %d \n",m,A->n);CHKERRQ(ierr); 329 ierr = PetscViewerASCIIPrintf(viewer,"%% Nonzeros = %d \n",a->nz);CHKERRQ(ierr); 330 ierr = PetscViewerASCIIPrintf(viewer,"zzz = zeros(%d,3);\n",a->nz+nofinalvalue);CHKERRQ(ierr); 331 ierr = PetscViewerASCIIPrintf(viewer,"zzz = [\n");CHKERRQ(ierr); 332 333 for (i=0; i<m; i++) { 334 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 335 #if defined(PETSC_USE_COMPLEX) 336 ierr = PetscViewerASCIIPrintf(viewer,"%d %d %18.16e + %18.16ei \n",i+1,a->j[j]+!shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 337 #else 338 ierr = PetscViewerASCIIPrintf(viewer,"%d %d %18.16e\n",i+1,a->j[j]+!shift,a->a[j]);CHKERRQ(ierr); 339 #endif 340 } 341 } 342 if (nofinalvalue) { 343 ierr = PetscViewerASCIIPrintf(viewer,"%d %d %18.16e\n",m,A->n,0.0);CHKERRQ(ierr); 344 } 345 ierr = PetscViewerASCIIPrintf(viewer,"];\n %s = spconvert(zzz);\n",name);CHKERRQ(ierr); 346 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 347 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 348 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 349 for (i=0; i<m; i++) { 350 ierr = PetscViewerASCIIPrintf(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 = PetscViewerASCIIPrintf(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 = PetscViewerASCIIPrintf(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 = PetscViewerASCIIPrintf(viewer," %d %g ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 359 } 360 #else 361 if (a->a[j] != 0.0) {ierr = PetscViewerASCIIPrintf(viewer," %d %g ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr);} 362 #endif 363 } 364 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 365 } 366 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 367 } else if (format == PETSC_VIEWER_ASCII_SYMMODU) { 368 int nzd=0,fshift=1,*sptr; 369 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 370 ierr = PetscMalloc((m+1)*sizeof(int),&sptr);CHKERRQ(ierr); 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 = PetscViewerASCIIPrintf(viewer," %d %d\n\n",m,nzd);CHKERRQ(ierr); 385 for (i=0; i<m+1; i+=6) { 386 if (i+4<m) {ierr = PetscViewerASCIIPrintf(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 = PetscViewerASCIIPrintf(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 = PetscViewerASCIIPrintf(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 = PetscViewerASCIIPrintf(viewer," %d %d %d\n",sptr[i],sptr[i+1],sptr[i+2]);CHKERRQ(ierr);} 390 else if (i<m) {ierr = PetscViewerASCIIPrintf(viewer," %d %d\n",sptr[i],sptr[i+1]);CHKERRQ(ierr);} 391 else {ierr = PetscViewerASCIIPrintf(viewer," %d\n",sptr[i]);CHKERRQ(ierr);} 392 } 393 ierr = PetscViewerASCIIPrintf(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 = PetscViewerASCIIPrintf(viewer," %d ",a->j[j]+fshift);CHKERRQ(ierr);} 398 } 399 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 400 } 401 ierr = PetscViewerASCIIPrintf(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 = PetscViewerASCIIPrintf(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 = PetscViewerASCIIPrintf(viewer," %18.16e ",a->a[j]);CHKERRQ(ierr);} 411 #endif 412 } 413 } 414 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 415 } 416 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 417 } else if (format == PETSC_VIEWER_ASCII_DENSE) { 418 int cnt = 0,jcnt; 419 Scalar value; 420 421 ierr = PetscViewerASCIIUseTabs(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 = PetscViewerASCIIPrintf(viewer," %7.5e+%7.5e i ",PetscRealPart(value),PetscImaginaryPart(value));CHKERRQ(ierr); 433 #else 434 ierr = PetscViewerASCIIPrintf(viewer," %7.5e ",value);CHKERRQ(ierr); 435 #endif 436 } 437 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 438 } 439 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 440 } else { 441 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 442 for (i=0; i<m; i++) { 443 ierr = PetscViewerASCIIPrintf(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 = PetscViewerASCIIPrintf(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 = PetscViewerASCIIPrintf(viewer," %d %g - %g i",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 450 } else { 451 ierr = PetscViewerASCIIPrintf(viewer," %d %g ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 452 } 453 #else 454 ierr = PetscViewerASCIIPrintf(viewer," %d %g ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr); 455 #endif 456 } 457 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 458 } 459 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 460 } 461 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 462 PetscFunctionReturn(0); 463 } 464 465 #undef __FUNC__ 466 #define __FUNC__ "MatView_SeqAIJ_Draw_Zoom" 467 int MatView_SeqAIJ_Draw_Zoom(PetscDraw 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 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r,maxv = 0.0; 473 PetscViewer viewer; 474 PetscViewerFormat format; 475 476 PetscFunctionBegin; 477 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 478 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 479 480 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 481 /* loop over matrix elements drawing boxes */ 482 483 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 484 /* Blue for negative, Cyan for zero and Red for positive */ 485 color = PETSC_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 = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 496 } 497 } 498 color = PETSC_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 = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 505 } 506 } 507 color = PETSC_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 = PetscDrawRectangle(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 PetscDraw 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 - PETSC_DRAW_BASIC_COLORS)/maxv; 531 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 532 if (popup) {ierr = PetscDrawScalePopup(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 = PETSC_DRAW_BASIC_COLORS + (int)(scale*PetscAbsScalar(a->a[count])); 539 ierr = PetscDrawRectangle(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__ "MatView_SeqAIJ_Draw" 549 int MatView_SeqAIJ_Draw(Mat A,PetscViewer viewer) 550 { 551 int ierr; 552 PetscDraw draw; 553 PetscReal xr,yr,xl,yl,h,w; 554 PetscTruth isnull; 555 556 PetscFunctionBegin; 557 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 558 ierr = PetscDrawIsNull(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 = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 565 ierr = PetscDrawZoom(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__ "MatView_SeqAIJ" 572 int MatView_SeqAIJ(Mat A,PetscViewer 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,PETSC_VIEWER_SOCKET,&issocket);CHKERRQ(ierr); 580 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);CHKERRQ(ierr); 581 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_BINARY,&isbinary);CHKERRQ(ierr); 582 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_DRAW,&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 = PetscViewerSocketPutSparse_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__ "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 PetscLogObjectMemory(A,-(m+1)*sizeof(int)); 643 a->diag = 0; 644 } 645 PetscLogInfo(A,"MatAssemblyEnd_SeqAIJ:Matrix size: %d X %d; storage space: %d unneeded,%d used\n",m,A->n,fshift,a->nz); 646 PetscLogInfo(A,"MatAssemblyEnd_SeqAIJ:Number of mallocs during MatSetValues() is %d\n",a->reallocs); 647 PetscLogInfo(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__ "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__ "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 PetscLogObjectState((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__ "MatCompress_SeqAIJ" 707 int MatCompress_SeqAIJ(Mat A) 708 { 709 PetscFunctionBegin; 710 PetscFunctionReturn(0); 711 } 712 713 #undef __FUNC__ 714 #define __FUNC__ "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_SYMMETRIC){ a->symmetric = PETSC_TRUE; 733 a->structurally_symmetric = PETSC_TRUE; 734 } 735 else if (op == MAT_STRUCTURALLY_SYMMETRIC) a->structurally_symmetric = PETSC_TRUE; 736 else if (op == MAT_ROWS_SORTED || 737 op == MAT_ROWS_UNSORTED || 738 op == MAT_YES_NEW_DIAGONALS || 739 op == MAT_IGNORE_OFF_PROC_ENTRIES|| 740 op == MAT_USE_HASH_TABLE) 741 PetscLogInfo(A,"MatSetOption_SeqAIJ:Option ignored\n"); 742 else if (op == MAT_NO_NEW_DIAGONALS) { 743 SETERRQ(PETSC_ERR_SUP,"MAT_NO_NEW_DIAGONALS"); 744 } else if (op == MAT_INODE_LIMIT_1) a->inode.limit = 1; 745 else if (op == MAT_INODE_LIMIT_2) a->inode.limit = 2; 746 else if (op == MAT_INODE_LIMIT_3) a->inode.limit = 3; 747 else if (op == MAT_INODE_LIMIT_4) a->inode.limit = 4; 748 else if (op == MAT_INODE_LIMIT_5) a->inode.limit = 5; 749 else SETERRQ(PETSC_ERR_SUP,"unknown option"); 750 PetscFunctionReturn(0); 751 } 752 753 #undef __FUNC__ 754 #define __FUNC__ "MatGetDiagonal_SeqAIJ" 755 int MatGetDiagonal_SeqAIJ(Mat A,Vec v) 756 { 757 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 758 int i,j,n,shift = a->indexshift,ierr; 759 Scalar *x,zero = 0.0; 760 761 PetscFunctionBegin; 762 ierr = VecSet(&zero,v);CHKERRQ(ierr); 763 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 764 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 765 if (n != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 766 for (i=0; i<A->m; i++) { 767 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 768 if (a->j[j]+shift == i) { 769 x[i] = a->a[j]; 770 break; 771 } 772 } 773 } 774 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 775 PetscFunctionReturn(0); 776 } 777 778 #undef __FUNC__ 779 #define __FUNC__ "MatMultTranspose_SeqAIJ" 780 int MatMultTranspose_SeqAIJ(Mat A,Vec xx,Vec yy) 781 { 782 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 783 Scalar *x,*y,*v,alpha,zero = 0.0; 784 int ierr,m = A->m,n,i,*idx,shift = a->indexshift; 785 786 PetscFunctionBegin; 787 ierr = VecSet(&zero,yy);CHKERRQ(ierr); 788 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 789 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 790 y = y + shift; /* shift for Fortran start by 1 indexing */ 791 for (i=0; i<m; i++) { 792 idx = a->j + a->i[i] + shift; 793 v = a->a + a->i[i] + shift; 794 n = a->i[i+1] - a->i[i]; 795 alpha = x[i]; 796 while (n-->0) {y[*idx++] += alpha * *v++;} 797 } 798 PetscLogFlops(2*a->nz - A->n); 799 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 800 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 801 PetscFunctionReturn(0); 802 } 803 804 #undef __FUNC__ 805 #define __FUNC__ "MatMultTransposeAdd_SeqAIJ" 806 int MatMultTransposeAdd_SeqAIJ(Mat A,Vec xx,Vec zz,Vec yy) 807 { 808 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 809 Scalar *x,*y,*v,alpha; 810 int ierr,m = A->m,n,i,*idx,shift = a->indexshift; 811 812 PetscFunctionBegin; 813 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 814 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 815 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 816 y = y + shift; /* shift for Fortran start by 1 indexing */ 817 for (i=0; i<m; i++) { 818 idx = a->j + a->i[i] + shift; 819 v = a->a + a->i[i] + shift; 820 n = a->i[i+1] - a->i[i]; 821 alpha = x[i]; 822 while (n-->0) {y[*idx++] += alpha * *v++;} 823 } 824 PetscLogFlops(2*a->nz); 825 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 826 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 827 PetscFunctionReturn(0); 828 } 829 830 #undef __FUNC__ 831 #define __FUNC__ "MatMult_SeqAIJ" 832 int MatMult_SeqAIJ(Mat A,Vec xx,Vec yy) 833 { 834 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 835 Scalar *x,*y,*v,sum; 836 int ierr,m = A->m,*idx,shift = a->indexshift,*ii; 837 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 838 int n,i,jrow,j; 839 #endif 840 841 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 842 #pragma disjoint(*x,*y,*v) 843 #endif 844 845 PetscFunctionBegin; 846 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 847 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 848 x = x + shift; /* shift for Fortran start by 1 indexing */ 849 idx = a->j; 850 v = a->a; 851 ii = a->i; 852 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 853 fortranmultaij_(&m,x,ii,idx+shift,v+shift,y); 854 #else 855 v += shift; /* shift for Fortran start by 1 indexing */ 856 idx += shift; 857 for (i=0; i<m; i++) { 858 jrow = ii[i]; 859 n = ii[i+1] - jrow; 860 sum = 0.0; 861 for (j=0; j<n; j++) { 862 sum += v[jrow]*x[idx[jrow]]; jrow++; 863 } 864 y[i] = sum; 865 } 866 #endif 867 PetscLogFlops(2*a->nz - m); 868 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 869 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 870 PetscFunctionReturn(0); 871 } 872 873 #undef __FUNC__ 874 #define __FUNC__ "MatMultAdd_SeqAIJ" 875 int MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 876 { 877 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 878 Scalar *x,*y,*z,*v,sum; 879 int ierr,m = A->m,*idx,shift = a->indexshift,*ii; 880 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 881 int n,i,jrow,j; 882 #endif 883 884 PetscFunctionBegin; 885 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 886 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 887 if (zz != yy) { 888 ierr = VecGetArray(zz,&z);CHKERRQ(ierr); 889 } else { 890 z = y; 891 } 892 x = x + shift; /* shift for Fortran start by 1 indexing */ 893 idx = a->j; 894 v = a->a; 895 ii = a->i; 896 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 897 fortranmultaddaij_(&m,x,ii,idx+shift,v+shift,y,z); 898 #else 899 v += shift; /* shift for Fortran start by 1 indexing */ 900 idx += shift; 901 for (i=0; i<m; i++) { 902 jrow = ii[i]; 903 n = ii[i+1] - jrow; 904 sum = y[i]; 905 for (j=0; j<n; j++) { 906 sum += v[jrow]*x[idx[jrow]]; jrow++; 907 } 908 z[i] = sum; 909 } 910 #endif 911 PetscLogFlops(2*a->nz); 912 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 913 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 914 if (zz != yy) { 915 ierr = VecRestoreArray(zz,&z);CHKERRQ(ierr); 916 } 917 PetscFunctionReturn(0); 918 } 919 920 /* 921 Adds diagonal pointers to sparse matrix structure. 922 */ 923 #undef __FUNC__ 924 #define __FUNC__ "MatMarkDiagonal_SeqAIJ" 925 int MatMarkDiagonal_SeqAIJ(Mat A) 926 { 927 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 928 int i,j,*diag,m = A->m,shift = a->indexshift,ierr; 929 930 PetscFunctionBegin; 931 if (a->diag) PetscFunctionReturn(0); 932 933 ierr = PetscMalloc((m+1)*sizeof(int),&diag);CHKERRQ(ierr); 934 PetscLogObjectMemory(A,(m+1)*sizeof(int)); 935 for (i=0; i<A->m; i++) { 936 diag[i] = a->i[i+1]; 937 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 938 if (a->j[j]+shift == i) { 939 diag[i] = j - shift; 940 break; 941 } 942 } 943 } 944 a->diag = diag; 945 PetscFunctionReturn(0); 946 } 947 948 /* 949 Checks for missing diagonals 950 */ 951 #undef __FUNC__ 952 #define __FUNC__ "MatMissingDiagonal_SeqAIJ" 953 int MatMissingDiagonal_SeqAIJ(Mat A) 954 { 955 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 956 int *diag,*jj = a->j,i,shift = a->indexshift,ierr; 957 958 PetscFunctionBegin; 959 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 960 diag = a->diag; 961 for (i=0; i<A->m; i++) { 962 if (jj[diag[i]+shift] != i-shift) { 963 SETERRQ1(1,"Matrix is missing diagonal number %d",i); 964 } 965 } 966 PetscFunctionReturn(0); 967 } 968 969 #undef __FUNC__ 970 #define __FUNC__ "MatRelax_SeqAIJ" 971 int MatRelax_SeqAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,int its,Vec xx) 972 { 973 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 974 Scalar *x,*b,*bs, d,*xs,sum,*v = a->a,*t=0,scale,*ts,*xb,*idiag=0; 975 int ierr,*idx,*diag,n = A->n,m = A->m,i,shift = a->indexshift; 976 977 PetscFunctionBegin; 978 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 979 if (xx != bb) { 980 ierr = VecGetArray(bb,&b);CHKERRQ(ierr); 981 } else { 982 b = x; 983 } 984 985 if (!a->diag) {ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr);} 986 diag = a->diag; 987 xs = x + shift; /* shifted by one for index start of a or a->j*/ 988 if (flag == SOR_APPLY_UPPER) { 989 /* apply (U + D/omega) to the vector */ 990 bs = b + shift; 991 for (i=0; i<m; i++) { 992 d = fshift + a->a[diag[i] + shift]; 993 n = a->i[i+1] - diag[i] - 1; 994 PetscLogFlops(2*n-1); 995 idx = a->j + diag[i] + (!shift); 996 v = a->a + diag[i] + (!shift); 997 sum = b[i]*d/omega; 998 SPARSEDENSEDOT(sum,bs,v,idx,n); 999 x[i] = sum; 1000 } 1001 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1002 if (bb != xx) {ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);} 1003 PetscFunctionReturn(0); 1004 } 1005 1006 /* setup workspace for Eisenstat */ 1007 if (flag & SOR_EISENSTAT) { 1008 if (!a->idiag) { 1009 ierr = PetscMalloc(2*m*sizeof(Scalar),&a->idiag);CHKERRQ(ierr); 1010 a->ssor = a->idiag + m; 1011 v = a->a; 1012 for (i=0; i<m; i++) { a->idiag[i] = 1.0/v[diag[i]];} 1013 } 1014 t = a->ssor; 1015 idiag = a->idiag; 1016 } 1017 /* Let A = L + U + D; where L is lower trianglar, 1018 U is upper triangular, E is diagonal; This routine applies 1019 1020 (L + E)^{-1} A (U + E)^{-1} 1021 1022 to a vector efficiently using Eisenstat's trick. This is for 1023 the case of SSOR preconditioner, so E is D/omega where omega 1024 is the relaxation factor. 1025 */ 1026 1027 if (flag == SOR_APPLY_LOWER) { 1028 SETERRQ(PETSC_ERR_SUP,"SOR_APPLY_LOWER is not implemented"); 1029 } else if ((flag & SOR_EISENSTAT) && omega == 1.0 && shift == 0 && fshift == 0.0) { 1030 /* special case for omega = 1.0 saves flops and some integer ops */ 1031 Scalar *v2; 1032 1033 v2 = a->a; 1034 /* x = (E + U)^{-1} b */ 1035 for (i=m-1; i>=0; i--) { 1036 n = a->i[i+1] - diag[i] - 1; 1037 idx = a->j + diag[i] + 1; 1038 v = a->a + diag[i] + 1; 1039 sum = b[i]; 1040 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1041 x[i] = sum*idiag[i]; 1042 1043 /* t = b - (2*E - D)x */ 1044 t[i] = b[i] - (v2[diag[i]])*x[i]; 1045 } 1046 1047 /* t = (E + L)^{-1}t */ 1048 diag = a->diag; 1049 for (i=0; i<m; i++) { 1050 n = diag[i] - a->i[i]; 1051 idx = a->j + a->i[i]; 1052 v = a->a + a->i[i]; 1053 sum = t[i]; 1054 SPARSEDENSEMDOT(sum,t,v,idx,n); 1055 t[i] = sum*idiag[i]; 1056 1057 /* x = x + t */ 1058 x[i] += t[i]; 1059 } 1060 1061 PetscLogFlops(3*m-1 + 2*a->nz); 1062 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1063 if (bb != xx) {ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);} 1064 PetscFunctionReturn(0); 1065 } else if (flag & SOR_EISENSTAT) { 1066 /* Let A = L + U + D; where L is lower trianglar, 1067 U is upper triangular, E is diagonal; This routine applies 1068 1069 (L + E)^{-1} A (U + E)^{-1} 1070 1071 to a vector efficiently using Eisenstat's trick. This is for 1072 the case of SSOR preconditioner, so E is D/omega where omega 1073 is the relaxation factor. 1074 */ 1075 scale = (2.0/omega) - 1.0; 1076 1077 /* x = (E + U)^{-1} b */ 1078 for (i=m-1; i>=0; i--) { 1079 d = fshift + a->a[diag[i] + shift]; 1080 n = a->i[i+1] - diag[i] - 1; 1081 idx = a->j + diag[i] + (!shift); 1082 v = a->a + diag[i] + (!shift); 1083 sum = b[i]; 1084 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1085 x[i] = omega*(sum/d); 1086 } 1087 1088 /* t = b - (2*E - D)x */ 1089 v = a->a; 1090 for (i=0; i<m; i++) { t[i] = b[i] - scale*(v[*diag++ + shift])*x[i]; } 1091 1092 /* t = (E + L)^{-1}t */ 1093 ts = t + shift; /* shifted by one for index start of a or a->j*/ 1094 diag = a->diag; 1095 for (i=0; i<m; i++) { 1096 d = fshift + a->a[diag[i]+shift]; 1097 n = diag[i] - a->i[i]; 1098 idx = a->j + a->i[i] + shift; 1099 v = a->a + a->i[i] + shift; 1100 sum = t[i]; 1101 SPARSEDENSEMDOT(sum,ts,v,idx,n); 1102 t[i] = omega*(sum/d); 1103 /* x = x + t */ 1104 x[i] += t[i]; 1105 } 1106 1107 PetscLogFlops(6*m-1 + 2*a->nz); 1108 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1109 if (bb != xx) {ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);} 1110 PetscFunctionReturn(0); 1111 } 1112 if (flag & SOR_ZERO_INITIAL_GUESS) { 1113 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1114 for (i=0; i<m; i++) { 1115 d = fshift + a->a[diag[i]+shift]; 1116 n = diag[i] - a->i[i]; 1117 PetscLogFlops(2*n-1); 1118 idx = a->j + a->i[i] + shift; 1119 v = a->a + a->i[i] + shift; 1120 sum = b[i]; 1121 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1122 x[i] = omega*(sum/d); 1123 } 1124 xb = x; 1125 } else xb = b; 1126 if ((flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) && 1127 (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP)) { 1128 for (i=0; i<m; i++) { 1129 x[i] *= a->a[diag[i]+shift]; 1130 } 1131 PetscLogFlops(m); 1132 } 1133 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1134 for (i=m-1; i>=0; i--) { 1135 d = fshift + a->a[diag[i] + shift]; 1136 n = a->i[i+1] - diag[i] - 1; 1137 PetscLogFlops(2*n-1); 1138 idx = a->j + diag[i] + (!shift); 1139 v = a->a + diag[i] + (!shift); 1140 sum = xb[i]; 1141 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1142 x[i] = omega*(sum/d); 1143 } 1144 } 1145 its--; 1146 } 1147 while (its--) { 1148 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1149 for (i=0; i<m; i++) { 1150 d = fshift + a->a[diag[i]+shift]; 1151 n = a->i[i+1] - a->i[i]; 1152 PetscLogFlops(2*n-1); 1153 idx = a->j + a->i[i] + shift; 1154 v = a->a + a->i[i] + shift; 1155 sum = b[i]; 1156 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1157 x[i] = (1. - omega)*x[i] + omega*(sum + a->a[diag[i]+shift]*x[i])/d; 1158 } 1159 } 1160 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1161 for (i=m-1; i>=0; i--) { 1162 d = fshift + a->a[diag[i] + shift]; 1163 n = a->i[i+1] - a->i[i]; 1164 PetscLogFlops(2*n-1); 1165 idx = a->j + a->i[i] + shift; 1166 v = a->a + a->i[i] + shift; 1167 sum = b[i]; 1168 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1169 x[i] = (1. - omega)*x[i] + omega*(sum + a->a[diag[i]+shift]*x[i])/d; 1170 } 1171 } 1172 } 1173 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1174 if (bb != xx) {ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);} 1175 PetscFunctionReturn(0); 1176 } 1177 1178 #undef __FUNC__ 1179 #define __FUNC__ "MatGetInfo_SeqAIJ" 1180 int MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1181 { 1182 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1183 1184 PetscFunctionBegin; 1185 info->rows_global = (double)A->m; 1186 info->columns_global = (double)A->n; 1187 info->rows_local = (double)A->m; 1188 info->columns_local = (double)A->n; 1189 info->block_size = 1.0; 1190 info->nz_allocated = (double)a->maxnz; 1191 info->nz_used = (double)a->nz; 1192 info->nz_unneeded = (double)(a->maxnz - a->nz); 1193 info->assemblies = (double)A->num_ass; 1194 info->mallocs = (double)a->reallocs; 1195 info->memory = A->mem; 1196 if (A->factor) { 1197 info->fill_ratio_given = A->info.fill_ratio_given; 1198 info->fill_ratio_needed = A->info.fill_ratio_needed; 1199 info->factor_mallocs = A->info.factor_mallocs; 1200 } else { 1201 info->fill_ratio_given = 0; 1202 info->fill_ratio_needed = 0; 1203 info->factor_mallocs = 0; 1204 } 1205 PetscFunctionReturn(0); 1206 } 1207 1208 EXTERN int MatLUFactorSymbolic_SeqAIJ(Mat,IS,IS,MatLUInfo*,Mat*); 1209 EXTERN int MatLUFactorNumeric_SeqAIJ(Mat,Mat*); 1210 EXTERN int MatLUFactor_SeqAIJ(Mat,IS,IS,MatLUInfo*); 1211 EXTERN int MatSolve_SeqAIJ(Mat,Vec,Vec); 1212 EXTERN int MatSolveAdd_SeqAIJ(Mat,Vec,Vec,Vec); 1213 EXTERN int MatSolveTranspose_SeqAIJ(Mat,Vec,Vec); 1214 EXTERN int MatSolveTransposeAdd_SeqAIJ(Mat,Vec,Vec,Vec); 1215 1216 #undef __FUNC__ 1217 #define __FUNC__ "MatZeroRows_SeqAIJ" 1218 int MatZeroRows_SeqAIJ(Mat A,IS is,Scalar *diag) 1219 { 1220 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1221 int i,ierr,N,*rows,m = A->m - 1,shift = a->indexshift; 1222 1223 PetscFunctionBegin; 1224 ierr = ISGetLocalSize(is,&N);CHKERRQ(ierr); 1225 ierr = ISGetIndices(is,&rows);CHKERRQ(ierr); 1226 if (a->keepzeroedrows) { 1227 for (i=0; i<N; i++) { 1228 if (rows[i] < 0 || rows[i] > m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"row out of range"); 1229 ierr = PetscMemzero(&a->a[a->i[rows[i]]+shift],a->ilen[rows[i]]*sizeof(Scalar));CHKERRQ(ierr); 1230 } 1231 if (diag) { 1232 ierr = MatMissingDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1233 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1234 for (i=0; i<N; i++) { 1235 a->a[a->diag[rows[i]]] = *diag; 1236 } 1237 } 1238 } else { 1239 if (diag) { 1240 for (i=0; i<N; i++) { 1241 if (rows[i] < 0 || rows[i] > m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"row out of range"); 1242 if (a->ilen[rows[i]] > 0) { 1243 a->ilen[rows[i]] = 1; 1244 a->a[a->i[rows[i]]+shift] = *diag; 1245 a->j[a->i[rows[i]]+shift] = rows[i]+shift; 1246 } else { /* in case row was completely empty */ 1247 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],diag,INSERT_VALUES);CHKERRQ(ierr); 1248 } 1249 } 1250 } else { 1251 for (i=0; i<N; i++) { 1252 if (rows[i] < 0 || rows[i] > m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"row out of range"); 1253 a->ilen[rows[i]] = 0; 1254 } 1255 } 1256 } 1257 ierr = ISRestoreIndices(is,&rows);CHKERRQ(ierr); 1258 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1259 PetscFunctionReturn(0); 1260 } 1261 1262 #undef __FUNC__ 1263 #define __FUNC__ "MatGetOwnershipRange_SeqAIJ" 1264 int MatGetOwnershipRange_SeqAIJ(Mat A,int *m,int *n) 1265 { 1266 PetscFunctionBegin; 1267 if (m) *m = 0; 1268 if (n) *n = A->m; 1269 PetscFunctionReturn(0); 1270 } 1271 1272 #undef __FUNC__ 1273 #define __FUNC__ "MatGetRow_SeqAIJ" 1274 int MatGetRow_SeqAIJ(Mat A,int row,int *nz,int **idx,Scalar **v) 1275 { 1276 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1277 int *itmp,i,shift = a->indexshift,ierr; 1278 1279 PetscFunctionBegin; 1280 if (row < 0 || row >= A->m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Row %d out of range",row); 1281 1282 *nz = a->i[row+1] - a->i[row]; 1283 if (v) *v = a->a + a->i[row] + shift; 1284 if (idx) { 1285 itmp = a->j + a->i[row] + shift; 1286 if (*nz && shift) { 1287 ierr = PetscMalloc((*nz)*sizeof(int),idx);CHKERRQ(ierr); 1288 for (i=0; i<(*nz); i++) {(*idx)[i] = itmp[i] + shift;} 1289 } else if (*nz) { 1290 *idx = itmp; 1291 } 1292 else *idx = 0; 1293 } 1294 PetscFunctionReturn(0); 1295 } 1296 1297 #undef __FUNC__ 1298 #define __FUNC__ "MatRestoreRow_SeqAIJ" 1299 int MatRestoreRow_SeqAIJ(Mat A,int row,int *nz,int **idx,Scalar **v) 1300 { 1301 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1302 int ierr; 1303 1304 PetscFunctionBegin; 1305 if (idx) {if (*idx && a->indexshift) {ierr = PetscFree(*idx);CHKERRQ(ierr);}} 1306 PetscFunctionReturn(0); 1307 } 1308 1309 #undef __FUNC__ 1310 #define __FUNC__ "MatNorm_SeqAIJ" 1311 int MatNorm_SeqAIJ(Mat A,NormType type,PetscReal *norm) 1312 { 1313 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1314 Scalar *v = a->a; 1315 PetscReal sum = 0.0; 1316 int i,j,shift = a->indexshift,ierr; 1317 1318 PetscFunctionBegin; 1319 if (type == NORM_FROBENIUS) { 1320 for (i=0; i<a->nz; i++) { 1321 #if defined(PETSC_USE_COMPLEX) 1322 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1323 #else 1324 sum += (*v)*(*v); v++; 1325 #endif 1326 } 1327 *norm = sqrt(sum); 1328 } else if (type == NORM_1) { 1329 PetscReal *tmp; 1330 int *jj = a->j; 1331 ierr = PetscMalloc((A->n+1)*sizeof(PetscReal),&tmp);CHKERRQ(ierr); 1332 ierr = PetscMemzero(tmp,A->n*sizeof(PetscReal));CHKERRQ(ierr); 1333 *norm = 0.0; 1334 for (j=0; j<a->nz; j++) { 1335 tmp[*jj++ + shift] += PetscAbsScalar(*v); v++; 1336 } 1337 for (j=0; j<A->n; j++) { 1338 if (tmp[j] > *norm) *norm = tmp[j]; 1339 } 1340 ierr = PetscFree(tmp);CHKERRQ(ierr); 1341 } else if (type == NORM_INFINITY) { 1342 *norm = 0.0; 1343 for (j=0; j<A->m; j++) { 1344 v = a->a + a->i[j] + shift; 1345 sum = 0.0; 1346 for (i=0; i<a->i[j+1]-a->i[j]; i++) { 1347 sum += PetscAbsScalar(*v); v++; 1348 } 1349 if (sum > *norm) *norm = sum; 1350 } 1351 } else { 1352 SETERRQ(PETSC_ERR_SUP,"No support for two norm"); 1353 } 1354 PetscFunctionReturn(0); 1355 } 1356 1357 #undef __FUNC__ 1358 #define __FUNC__ "MatTranspose_SeqAIJ" 1359 int MatTranspose_SeqAIJ(Mat A,Mat *B) 1360 { 1361 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1362 Mat C; 1363 int i,ierr,*aj = a->j,*ai = a->i,m = A->m,len,*col; 1364 int shift = a->indexshift; 1365 Scalar *array = a->a; 1366 1367 PetscFunctionBegin; 1368 if (!B && m != A->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Square matrix only for in-place"); 1369 ierr = PetscMalloc((1+A->n)*sizeof(int),&col);CHKERRQ(ierr); 1370 ierr = PetscMemzero(col,(1+A->n)*sizeof(int));CHKERRQ(ierr); 1371 if (shift) { 1372 for (i=0; i<ai[m]-1; i++) aj[i] -= 1; 1373 } 1374 for (i=0; i<ai[m]+shift; i++) col[aj[i]] += 1; 1375 ierr = MatCreateSeqAIJ(A->comm,A->n,m,0,col,&C);CHKERRQ(ierr); 1376 ierr = PetscFree(col);CHKERRQ(ierr); 1377 for (i=0; i<m; i++) { 1378 len = ai[i+1]-ai[i]; 1379 ierr = MatSetValues(C,len,aj,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 1380 array += len; 1381 aj += len; 1382 } 1383 if (shift) { 1384 for (i=0; i<ai[m]-1; i++) aj[i] += 1; 1385 } 1386 1387 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1388 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1389 1390 if (B) { 1391 *B = C; 1392 } else { 1393 ierr = MatHeaderCopy(A,C);CHKERRQ(ierr); 1394 } 1395 PetscFunctionReturn(0); 1396 } 1397 1398 #undef __FUNC__ 1399 #define __FUNC__ "MatDiagonalScale_SeqAIJ" 1400 int MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 1401 { 1402 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1403 Scalar *l,*r,x,*v; 1404 int ierr,i,j,m = A->m,n = A->n,M,nz = a->nz,*jj,shift = a->indexshift; 1405 1406 PetscFunctionBegin; 1407 if (ll) { 1408 /* The local size is used so that VecMPI can be passed to this routine 1409 by MatDiagonalScale_MPIAIJ */ 1410 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 1411 if (m != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 1412 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1413 v = a->a; 1414 for (i=0; i<m; i++) { 1415 x = l[i]; 1416 M = a->i[i+1] - a->i[i]; 1417 for (j=0; j<M; j++) { (*v++) *= x;} 1418 } 1419 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1420 PetscLogFlops(nz); 1421 } 1422 if (rr) { 1423 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 1424 if (n != A->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 1425 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1426 v = a->a; jj = a->j; 1427 for (i=0; i<nz; i++) { 1428 (*v++) *= r[*jj++ + shift]; 1429 } 1430 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1431 PetscLogFlops(nz); 1432 } 1433 PetscFunctionReturn(0); 1434 } 1435 1436 #undef __FUNC__ 1437 #define __FUNC__ "MatGetSubMatrix_SeqAIJ" 1438 int MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,int csize,MatReuse scall,Mat *B) 1439 { 1440 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 1441 int *smap,i,k,kstart,kend,ierr,oldcols = A->n,*lens; 1442 int row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 1443 int *irow,*icol,nrows,ncols,shift = a->indexshift,*ssmap; 1444 int *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 1445 Scalar *a_new,*mat_a; 1446 Mat C; 1447 PetscTruth stride; 1448 1449 PetscFunctionBegin; 1450 ierr = ISSorted(isrow,(PetscTruth*)&i);CHKERRQ(ierr); 1451 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"ISrow is not sorted"); 1452 ierr = ISSorted(iscol,(PetscTruth*)&i);CHKERRQ(ierr); 1453 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"IScol is not sorted"); 1454 1455 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1456 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1457 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1458 1459 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 1460 ierr = ISStride(iscol,&stride);CHKERRQ(ierr); 1461 if (stride && step == 1) { 1462 /* special case of contiguous rows */ 1463 ierr = PetscMalloc((ncols+nrows+1)*sizeof(int),&lens);CHKERRQ(ierr); 1464 starts = lens + ncols; 1465 /* loop over new rows determining lens and starting points */ 1466 for (i=0; i<nrows; i++) { 1467 kstart = ai[irow[i]]+shift; 1468 kend = kstart + ailen[irow[i]]; 1469 for (k=kstart; k<kend; k++) { 1470 if (aj[k]+shift >= first) { 1471 starts[i] = k; 1472 break; 1473 } 1474 } 1475 sum = 0; 1476 while (k < kend) { 1477 if (aj[k++]+shift >= first+ncols) break; 1478 sum++; 1479 } 1480 lens[i] = sum; 1481 } 1482 /* create submatrix */ 1483 if (scall == MAT_REUSE_MATRIX) { 1484 int n_cols,n_rows; 1485 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1486 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 1487 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 1488 C = *B; 1489 } else { 1490 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1491 } 1492 c = (Mat_SeqAIJ*)C->data; 1493 1494 /* loop over rows inserting into submatrix */ 1495 a_new = c->a; 1496 j_new = c->j; 1497 i_new = c->i; 1498 i_new[0] = -shift; 1499 for (i=0; i<nrows; i++) { 1500 ii = starts[i]; 1501 lensi = lens[i]; 1502 for (k=0; k<lensi; k++) { 1503 *j_new++ = aj[ii+k] - first; 1504 } 1505 ierr = PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(Scalar));CHKERRQ(ierr); 1506 a_new += lensi; 1507 i_new[i+1] = i_new[i] + lensi; 1508 c->ilen[i] = lensi; 1509 } 1510 ierr = PetscFree(lens);CHKERRQ(ierr); 1511 } else { 1512 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1513 ierr = PetscMalloc((1+oldcols)*sizeof(int),&smap);CHKERRQ(ierr); 1514 ssmap = smap + shift; 1515 ierr = PetscMalloc((1+nrows)*sizeof(int),&lens);CHKERRQ(ierr); 1516 ierr = PetscMemzero(smap,oldcols*sizeof(int));CHKERRQ(ierr); 1517 for (i=0; i<ncols; i++) smap[icol[i]] = i+1; 1518 /* determine lens of each row */ 1519 for (i=0; i<nrows; i++) { 1520 kstart = ai[irow[i]]+shift; 1521 kend = kstart + a->ilen[irow[i]]; 1522 lens[i] = 0; 1523 for (k=kstart; k<kend; k++) { 1524 if (ssmap[aj[k]]) { 1525 lens[i]++; 1526 } 1527 } 1528 } 1529 /* Create and fill new matrix */ 1530 if (scall == MAT_REUSE_MATRIX) { 1531 PetscTruth equal; 1532 1533 c = (Mat_SeqAIJ *)((*B)->data); 1534 if ((*B)->m != nrows || (*B)->n != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 1535 ierr = PetscMemcmp(c->ilen,lens,(*B)->m*sizeof(int),&equal);CHKERRQ(ierr); 1536 if (!equal) { 1537 SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 1538 } 1539 ierr = PetscMemzero(c->ilen,(*B)->m*sizeof(int));CHKERRQ(ierr); 1540 C = *B; 1541 } else { 1542 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1543 } 1544 c = (Mat_SeqAIJ *)(C->data); 1545 for (i=0; i<nrows; i++) { 1546 row = irow[i]; 1547 kstart = ai[row]+shift; 1548 kend = kstart + a->ilen[row]; 1549 mat_i = c->i[i]+shift; 1550 mat_j = c->j + mat_i; 1551 mat_a = c->a + mat_i; 1552 mat_ilen = c->ilen + i; 1553 for (k=kstart; k<kend; k++) { 1554 if ((tcol=ssmap[a->j[k]])) { 1555 *mat_j++ = tcol - (!shift); 1556 *mat_a++ = a->a[k]; 1557 (*mat_ilen)++; 1558 1559 } 1560 } 1561 } 1562 /* Free work space */ 1563 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1564 ierr = PetscFree(smap);CHKERRQ(ierr); 1565 ierr = PetscFree(lens);CHKERRQ(ierr); 1566 } 1567 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1568 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1569 1570 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1571 *B = C; 1572 PetscFunctionReturn(0); 1573 } 1574 1575 /* 1576 */ 1577 #undef __FUNC__ 1578 #define __FUNC__ "MatILUFactor_SeqAIJ" 1579 int MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,MatILUInfo *info) 1580 { 1581 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1582 int ierr; 1583 Mat outA; 1584 PetscTruth row_identity,col_identity; 1585 1586 PetscFunctionBegin; 1587 if (info && info->levels != 0) SETERRQ(PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 1588 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 1589 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 1590 if (!row_identity || !col_identity) { 1591 SETERRQ(1,"Row and column permutations must be identity for in-place ILU"); 1592 } 1593 1594 outA = inA; 1595 inA->factor = FACTOR_LU; 1596 a->row = row; 1597 a->col = col; 1598 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1599 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 1600 1601 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 1602 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} /* need to remove old one */ 1603 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 1604 PetscLogObjectParent(inA,a->icol); 1605 1606 if (!a->solve_work) { /* this matrix may have been factored before */ 1607 ierr = PetscMalloc((inA->m+1)*sizeof(Scalar),&a->solve_work);CHKERRQ(ierr); 1608 } 1609 1610 if (!a->diag) { 1611 ierr = MatMarkDiagonal_SeqAIJ(inA);CHKERRQ(ierr); 1612 } 1613 ierr = MatLUFactorNumeric_SeqAIJ(inA,&outA);CHKERRQ(ierr); 1614 PetscFunctionReturn(0); 1615 } 1616 1617 #include "petscblaslapack.h" 1618 #undef __FUNC__ 1619 #define __FUNC__ "MatScale_SeqAIJ" 1620 int MatScale_SeqAIJ(Scalar *alpha,Mat inA) 1621 { 1622 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1623 int one = 1; 1624 1625 PetscFunctionBegin; 1626 BLscal_(&a->nz,alpha,a->a,&one); 1627 PetscLogFlops(a->nz); 1628 PetscFunctionReturn(0); 1629 } 1630 1631 #undef __FUNC__ 1632 #define __FUNC__ "MatGetSubMatrices_SeqAIJ" 1633 int MatGetSubMatrices_SeqAIJ(Mat A,int n,IS *irow,IS *icol,MatReuse scall,Mat **B) 1634 { 1635 int ierr,i; 1636 1637 PetscFunctionBegin; 1638 if (scall == MAT_INITIAL_MATRIX) { 1639 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1640 } 1641 1642 for (i=0; i<n; i++) { 1643 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1644 } 1645 PetscFunctionReturn(0); 1646 } 1647 1648 #undef __FUNC__ 1649 #define __FUNC__ "MatGetBlockSize_SeqAIJ" 1650 int MatGetBlockSize_SeqAIJ(Mat A,int *bs) 1651 { 1652 PetscFunctionBegin; 1653 *bs = 1; 1654 PetscFunctionReturn(0); 1655 } 1656 1657 #undef __FUNC__ 1658 #define __FUNC__ "MatIncreaseOverlap_SeqAIJ" 1659 int MatIncreaseOverlap_SeqAIJ(Mat A,int is_max,IS *is,int ov) 1660 { 1661 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1662 int shift,row,i,j,k,l,m,n,*idx,ierr,*nidx,isz,val; 1663 int start,end,*ai,*aj; 1664 PetscBT table; 1665 1666 PetscFunctionBegin; 1667 shift = a->indexshift; 1668 m = A->m; 1669 ai = a->i; 1670 aj = a->j+shift; 1671 1672 if (ov < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"illegal overlap value used"); 1673 1674 ierr = PetscMalloc((m+1)*sizeof(int),&nidx);CHKERRQ(ierr); 1675 ierr = PetscBTCreate(m,table);CHKERRQ(ierr); 1676 1677 for (i=0; i<is_max; i++) { 1678 /* Initialize the two local arrays */ 1679 isz = 0; 1680 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 1681 1682 /* Extract the indices, assume there can be duplicate entries */ 1683 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 1684 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 1685 1686 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 1687 for (j=0; j<n ; ++j){ 1688 if(!PetscBTLookupSet(table,idx[j])) { nidx[isz++] = idx[j];} 1689 } 1690 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 1691 ierr = ISDestroy(is[i]);CHKERRQ(ierr); 1692 1693 k = 0; 1694 for (j=0; j<ov; j++){ /* for each overlap */ 1695 n = isz; 1696 for (; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 1697 row = nidx[k]; 1698 start = ai[row]; 1699 end = ai[row+1]; 1700 for (l = start; l<end ; l++){ 1701 val = aj[l] + shift; 1702 if (!PetscBTLookupSet(table,val)) {nidx[isz++] = val;} 1703 } 1704 } 1705 } 1706 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,(is+i));CHKERRQ(ierr); 1707 } 1708 ierr = PetscBTDestroy(table);CHKERRQ(ierr); 1709 ierr = PetscFree(nidx);CHKERRQ(ierr); 1710 PetscFunctionReturn(0); 1711 } 1712 1713 /* -------------------------------------------------------------- */ 1714 #undef __FUNC__ 1715 #define __FUNC__ "MatPermute_SeqAIJ" 1716 int MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 1717 { 1718 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1719 Scalar *vwork; 1720 int i,ierr,nz,m = A->m,n = A->n,*cwork; 1721 int *row,*col,*cnew,j,*lens; 1722 IS icolp,irowp; 1723 1724 PetscFunctionBegin; 1725 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 1726 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 1727 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 1728 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 1729 1730 /* determine lengths of permuted rows */ 1731 ierr = PetscMalloc((m+1)*sizeof(int),&lens);CHKERRQ(ierr); 1732 for (i=0; i<m; i++) { 1733 lens[row[i]] = a->i[i+1] - a->i[i]; 1734 } 1735 ierr = MatCreateSeqAIJ(A->comm,m,n,0,lens,B);CHKERRQ(ierr); 1736 ierr = PetscFree(lens);CHKERRQ(ierr); 1737 1738 ierr = PetscMalloc(n*sizeof(int),&cnew);CHKERRQ(ierr); 1739 for (i=0; i<m; i++) { 1740 ierr = MatGetRow(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 1741 for (j=0; j<nz; j++) { cnew[j] = col[cwork[j]];} 1742 ierr = MatSetValues(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 1743 ierr = MatRestoreRow(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 1744 } 1745 ierr = PetscFree(cnew);CHKERRQ(ierr); 1746 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1747 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1748 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 1749 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 1750 ierr = ISDestroy(irowp);CHKERRQ(ierr); 1751 ierr = ISDestroy(icolp);CHKERRQ(ierr); 1752 PetscFunctionReturn(0); 1753 } 1754 1755 #undef __FUNC__ 1756 #define __FUNC__ "MatPrintHelp_SeqAIJ" 1757 int MatPrintHelp_SeqAIJ(Mat A) 1758 { 1759 static PetscTruth called = PETSC_FALSE; 1760 MPI_Comm comm = A->comm; 1761 int ierr; 1762 1763 PetscFunctionBegin; 1764 if (called) {PetscFunctionReturn(0);} else called = PETSC_TRUE; 1765 ierr = (*PetscHelpPrintf)(comm," Options for MATSEQAIJ and MATMPIAIJ matrix formats (the defaults):\n");CHKERRQ(ierr); 1766 ierr = (*PetscHelpPrintf)(comm," -mat_lu_pivotthreshold <threshold>: Set pivoting threshold\n");CHKERRQ(ierr); 1767 ierr = (*PetscHelpPrintf)(comm," -mat_aij_oneindex: internal indices begin at 1 instead of the default 0.\n");CHKERRQ(ierr); 1768 ierr = (*PetscHelpPrintf)(comm," -mat_aij_no_inode: Do not use inodes\n");CHKERRQ(ierr); 1769 ierr = (*PetscHelpPrintf)(comm," -mat_aij_inode_limit <limit>: Set inode limit (max limit=5)\n");CHKERRQ(ierr); 1770 #if defined(PETSC_HAVE_ESSL) 1771 ierr = (*PetscHelpPrintf)(comm," -mat_aij_essl: Use IBM sparse LU factorization and solve.\n");CHKERRQ(ierr); 1772 #endif 1773 #if defined(PETSC_HAVE_LUSOL) 1774 ierr = (*PetscHelpPrintf)(comm," -mat_aij_lusol: Use the Stanford LUSOL sparse factorization and solve.\n");CHKERRQ(ierr); 1775 #endif 1776 #if defined(PETSC_HAVE_MATLAB_ENGINE) && !defined(PETSC_USE_COMPLEX) 1777 ierr = (*PetscHelpPrintf)(comm," -mat_aij_matlab: Use Matlab engine sparse LU factorization and solve.\n");CHKERRQ(ierr); 1778 #endif 1779 PetscFunctionReturn(0); 1780 } 1781 EXTERN int MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg); 1782 EXTERN int MatFDColoringCreate_SeqAIJ(Mat,ISColoring,MatFDColoring); 1783 EXTERN int MatColoringPatch_SeqAIJ(Mat,int,int *,ISColoring *); 1784 EXTERN int MatILUDTFactor_SeqAIJ(Mat,MatILUInfo*,IS,IS,Mat*); 1785 #undef __FUNC__ 1786 #define __FUNC__ "MatCopy_SeqAIJ" 1787 int MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 1788 { 1789 int ierr; 1790 PetscTruth flg; 1791 1792 PetscFunctionBegin; 1793 ierr = PetscTypeCompare((PetscObject)B,MATSEQAIJ,&flg);CHKERRQ(ierr); 1794 if (str == SAME_NONZERO_PATTERN && flg) { 1795 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1796 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 1797 1798 if (a->i[A->m]+a->indexshift != b->i[B->m]+a->indexshift) { 1799 SETERRQ(1,"Number of nonzeros in two matrices are different"); 1800 } 1801 ierr = PetscMemcpy(b->a,a->a,(a->i[A->m]+a->indexshift)*sizeof(Scalar));CHKERRQ(ierr); 1802 } else { 1803 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1804 } 1805 PetscFunctionReturn(0); 1806 } 1807 1808 #undef __FUNC__ 1809 #define __FUNC__ "MatSetUpPreallocation_SeqAIJ" 1810 int MatSetUpPreallocation_SeqAIJ(Mat A) 1811 { 1812 int ierr; 1813 1814 PetscFunctionBegin; 1815 ierr = MatSeqAIJSetPreallocation(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 1816 PetscFunctionReturn(0); 1817 } 1818 1819 #undef __FUNC__ 1820 #define __FUNC__ "MatGetArray_SeqAIJ" 1821 int MatGetArray_SeqAIJ(Mat A,Scalar **array) 1822 { 1823 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1824 PetscFunctionBegin; 1825 *array = a->a; 1826 PetscFunctionReturn(0); 1827 } 1828 1829 #undef __FUNC__ 1830 #define __FUNC__ "MatRestoreArray_SeqAIJ" 1831 int MatRestoreArray_SeqAIJ(Mat A,Scalar **array) 1832 { 1833 PetscFunctionBegin; 1834 PetscFunctionReturn(0); 1835 } 1836 1837 /* -------------------------------------------------------------------*/ 1838 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 1839 MatGetRow_SeqAIJ, 1840 MatRestoreRow_SeqAIJ, 1841 MatMult_SeqAIJ, 1842 MatMultAdd_SeqAIJ, 1843 MatMultTranspose_SeqAIJ, 1844 MatMultTransposeAdd_SeqAIJ, 1845 MatSolve_SeqAIJ, 1846 MatSolveAdd_SeqAIJ, 1847 MatSolveTranspose_SeqAIJ, 1848 MatSolveTransposeAdd_SeqAIJ, 1849 MatLUFactor_SeqAIJ, 1850 0, 1851 MatRelax_SeqAIJ, 1852 MatTranspose_SeqAIJ, 1853 MatGetInfo_SeqAIJ, 1854 MatEqual_SeqAIJ, 1855 MatGetDiagonal_SeqAIJ, 1856 MatDiagonalScale_SeqAIJ, 1857 MatNorm_SeqAIJ, 1858 0, 1859 MatAssemblyEnd_SeqAIJ, 1860 MatCompress_SeqAIJ, 1861 MatSetOption_SeqAIJ, 1862 MatZeroEntries_SeqAIJ, 1863 MatZeroRows_SeqAIJ, 1864 MatLUFactorSymbolic_SeqAIJ, 1865 MatLUFactorNumeric_SeqAIJ, 1866 0, 1867 0, 1868 MatSetUpPreallocation_SeqAIJ, 1869 0, 1870 MatGetOwnershipRange_SeqAIJ, 1871 MatILUFactorSymbolic_SeqAIJ, 1872 0, 1873 MatGetArray_SeqAIJ, 1874 MatRestoreArray_SeqAIJ, 1875 MatDuplicate_SeqAIJ, 1876 0, 1877 0, 1878 MatILUFactor_SeqAIJ, 1879 0, 1880 0, 1881 MatGetSubMatrices_SeqAIJ, 1882 MatIncreaseOverlap_SeqAIJ, 1883 MatGetValues_SeqAIJ, 1884 MatCopy_SeqAIJ, 1885 MatPrintHelp_SeqAIJ, 1886 MatScale_SeqAIJ, 1887 0, 1888 0, 1889 MatILUDTFactor_SeqAIJ, 1890 MatGetBlockSize_SeqAIJ, 1891 MatGetRowIJ_SeqAIJ, 1892 MatRestoreRowIJ_SeqAIJ, 1893 MatGetColumnIJ_SeqAIJ, 1894 MatRestoreColumnIJ_SeqAIJ, 1895 MatFDColoringCreate_SeqAIJ, 1896 MatColoringPatch_SeqAIJ, 1897 0, 1898 MatPermute_SeqAIJ, 1899 0, 1900 0, 1901 MatDestroy_SeqAIJ, 1902 MatView_SeqAIJ, 1903 MatGetMaps_Petsc}; 1904 1905 EXTERN int MatUseSuperLU_SeqAIJ(Mat); 1906 EXTERN int MatUseEssl_SeqAIJ(Mat); 1907 EXTERN int MatUseLUSOL_SeqAIJ(Mat); 1908 EXTERN int MatUseMatlab_SeqAIJ(Mat); 1909 EXTERN int MatUseDXML_SeqAIJ(Mat); 1910 1911 EXTERN_C_BEGIN 1912 #undef __FUNC__ 1913 #define __FUNC__ "MatSeqAIJSetColumnIndices_SeqAIJ" 1914 1915 int MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,int *indices) 1916 { 1917 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 1918 int i,nz,n; 1919 1920 PetscFunctionBegin; 1921 if (aij->indexshift) SETERRQ(1,"No support with 1 based indexing"); 1922 1923 nz = aij->maxnz; 1924 n = mat->n; 1925 for (i=0; i<nz; i++) { 1926 aij->j[i] = indices[i]; 1927 } 1928 aij->nz = nz; 1929 for (i=0; i<n; i++) { 1930 aij->ilen[i] = aij->imax[i]; 1931 } 1932 1933 PetscFunctionReturn(0); 1934 } 1935 EXTERN_C_END 1936 1937 #undef __FUNC__ 1938 #define __FUNC__ "MatSeqAIJSetColumnIndices" 1939 /*@ 1940 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 1941 in the matrix. 1942 1943 Input Parameters: 1944 + mat - the SeqAIJ matrix 1945 - indices - the column indices 1946 1947 Level: advanced 1948 1949 Notes: 1950 This can be called if you have precomputed the nonzero structure of the 1951 matrix and want to provide it to the matrix object to improve the performance 1952 of the MatSetValues() operation. 1953 1954 You MUST have set the correct numbers of nonzeros per row in the call to 1955 MatCreateSeqAIJ(). 1956 1957 MUST be called before any calls to MatSetValues(); 1958 1959 @*/ 1960 int MatSeqAIJSetColumnIndices(Mat mat,int *indices) 1961 { 1962 int ierr,(*f)(Mat,int *); 1963 1964 PetscFunctionBegin; 1965 PetscValidHeaderSpecific(mat,MAT_COOKIE); 1966 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqAIJSetColumnIndices_C",(void **)&f);CHKERRQ(ierr); 1967 if (f) { 1968 ierr = (*f)(mat,indices);CHKERRQ(ierr); 1969 } else { 1970 SETERRQ(1,"Wrong type of matrix to set column indices"); 1971 } 1972 PetscFunctionReturn(0); 1973 } 1974 1975 /* ----------------------------------------------------------------------------------------*/ 1976 1977 EXTERN_C_BEGIN 1978 #undef __FUNC__ 1979 #define __FUNC__ "MatStoreValues_SeqAIJ" 1980 int MatStoreValues_SeqAIJ(Mat mat) 1981 { 1982 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 1983 int nz = aij->i[mat->m]+aij->indexshift,ierr; 1984 1985 PetscFunctionBegin; 1986 if (aij->nonew != 1) { 1987 SETERRQ(1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 1988 } 1989 1990 /* allocate space for values if not already there */ 1991 if (!aij->saved_values) { 1992 ierr = PetscMalloc((nz+1)*sizeof(Scalar),&aij->saved_values);CHKERRQ(ierr); 1993 } 1994 1995 /* copy values over */ 1996 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(Scalar));CHKERRQ(ierr); 1997 PetscFunctionReturn(0); 1998 } 1999 EXTERN_C_END 2000 2001 #undef __FUNC__ 2002 #define __FUNC__ /*<a name="MatStoreValues""></a>*/"MatStoreValues" 2003 /*@ 2004 MatStoreValues - Stashes a copy of the matrix values; this allows, for 2005 example, reuse of the linear part of a Jacobian, while recomputing the 2006 nonlinear portion. 2007 2008 Collect on Mat 2009 2010 Input Parameters: 2011 . mat - the matrix (currently on AIJ matrices support this option) 2012 2013 Level: advanced 2014 2015 Common Usage, with SNESSolve(): 2016 $ Create Jacobian matrix 2017 $ Set linear terms into matrix 2018 $ Apply boundary conditions to matrix, at this time matrix must have 2019 $ final nonzero structure (i.e. setting the nonlinear terms and applying 2020 $ boundary conditions again will not change the nonzero structure 2021 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 2022 $ ierr = MatStoreValues(mat); 2023 $ Call SNESSetJacobian() with matrix 2024 $ In your Jacobian routine 2025 $ ierr = MatRetrieveValues(mat); 2026 $ Set nonlinear terms in matrix 2027 2028 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 2029 $ // build linear portion of Jacobian 2030 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 2031 $ ierr = MatStoreValues(mat); 2032 $ loop over nonlinear iterations 2033 $ ierr = MatRetrieveValues(mat); 2034 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 2035 $ // call MatAssemblyBegin/End() on matrix 2036 $ Solve linear system with Jacobian 2037 $ endloop 2038 2039 Notes: 2040 Matrix must already be assemblied before calling this routine 2041 Must set the matrix option MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); before 2042 calling this routine. 2043 2044 .seealso: MatRetrieveValues() 2045 2046 @*/ 2047 int MatStoreValues(Mat mat) 2048 { 2049 int ierr,(*f)(Mat); 2050 2051 PetscFunctionBegin; 2052 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2053 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2054 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2055 2056 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatStoreValues_C",(void **)&f);CHKERRQ(ierr); 2057 if (f) { 2058 ierr = (*f)(mat);CHKERRQ(ierr); 2059 } else { 2060 SETERRQ(1,"Wrong type of matrix to store values"); 2061 } 2062 PetscFunctionReturn(0); 2063 } 2064 2065 EXTERN_C_BEGIN 2066 #undef __FUNC__ 2067 #define __FUNC__ "MatRetrieveValues_SeqAIJ" 2068 int MatRetrieveValues_SeqAIJ(Mat mat) 2069 { 2070 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2071 int nz = aij->i[mat->m]+aij->indexshift,ierr; 2072 2073 PetscFunctionBegin; 2074 if (aij->nonew != 1) { 2075 SETERRQ(1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 2076 } 2077 if (!aij->saved_values) { 2078 SETERRQ(1,"Must call MatStoreValues(A);first"); 2079 } 2080 2081 /* copy values over */ 2082 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(Scalar));CHKERRQ(ierr); 2083 PetscFunctionReturn(0); 2084 } 2085 EXTERN_C_END 2086 2087 #undef __FUNC__ 2088 #define __FUNC__ "MatRetrieveValues" 2089 /*@ 2090 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 2091 example, reuse of the linear part of a Jacobian, while recomputing the 2092 nonlinear portion. 2093 2094 Collect on Mat 2095 2096 Input Parameters: 2097 . mat - the matrix (currently on AIJ matrices support this option) 2098 2099 Level: advanced 2100 2101 .seealso: MatStoreValues() 2102 2103 @*/ 2104 int MatRetrieveValues(Mat mat) 2105 { 2106 int ierr,(*f)(Mat); 2107 2108 PetscFunctionBegin; 2109 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2110 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2111 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2112 2113 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatRetrieveValues_C",(void **)&f);CHKERRQ(ierr); 2114 if (f) { 2115 ierr = (*f)(mat);CHKERRQ(ierr); 2116 } else { 2117 SETERRQ(1,"Wrong type of matrix to retrieve values"); 2118 } 2119 PetscFunctionReturn(0); 2120 } 2121 2122 /* 2123 This allows SeqAIJ matrices to be passed to the matlab engine 2124 */ 2125 #if defined(PETSC_HAVE_MATLAB_ENGINE) && !defined(PETSC_USE_COMPLEX) 2126 #include "engine.h" /* Matlab include file */ 2127 #include "mex.h" /* Matlab include file */ 2128 EXTERN_C_BEGIN 2129 #undef __FUNC__ 2130 #define __FUNC__ "MatMatlabEnginePut_SeqAIJ" 2131 int MatMatlabEnginePut_SeqAIJ(PetscObject obj,void *engine) 2132 { 2133 int ierr,i,*ai,*aj; 2134 Mat B = (Mat)obj; 2135 Scalar *array; 2136 mxArray *mat; 2137 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)B->data; 2138 2139 PetscFunctionBegin; 2140 mat = mxCreateSparse(B->n,B->m,aij->nz,mxREAL); 2141 ierr = PetscMemcpy(mxGetPr(mat),aij->a,aij->nz*sizeof(Scalar));CHKERRQ(ierr); 2142 /* Matlab stores by column, not row so we pass in the transpose of the matrix */ 2143 ierr = PetscMemcpy(mxGetIr(mat),aij->j,aij->nz*sizeof(int));CHKERRQ(ierr); 2144 ierr = PetscMemcpy(mxGetJc(mat),aij->i,(B->m+1)*sizeof(int));CHKERRQ(ierr); 2145 2146 /* Matlab indices start at 0 for sparse (what a surprise) */ 2147 if (aij->indexshift) { 2148 for (i=0; i<B->m+1; i++) { 2149 ai[i]--; 2150 } 2151 for (i=0; i<aij->nz; i++) { 2152 aj[i]--; 2153 } 2154 } 2155 ierr = PetscObjectName(obj);CHKERRQ(ierr); 2156 mxSetName(mat,obj->name); 2157 engPutArray((Engine *)engine,mat); 2158 PetscFunctionReturn(0); 2159 } 2160 EXTERN_C_END 2161 2162 EXTERN_C_BEGIN 2163 #undef __FUNC__ 2164 #define __FUNC__ "MatMatlabEngineGet_SeqAIJ" 2165 int MatMatlabEngineGet_SeqAIJ(PetscObject obj,void *engine) 2166 { 2167 int ierr,ii; 2168 Mat mat = (Mat)obj; 2169 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 2170 mxArray *mmat; 2171 2172 PetscFunctionBegin; 2173 ierr = PetscFree(aij->a);CHKERRQ(ierr); 2174 aij->indexshift = 0; 2175 2176 mmat = engGetArray((Engine *)engine,obj->name); 2177 2178 aij->nz = (mxGetJc(mmat))[mat->m]; 2179 ierr = PetscMalloc(aij->nz*(sizeof(int)+sizeof(Scalar))+(mat->m+1)*sizeof(int),&aij->a);CHKERRQ(ierr); 2180 aij->j = (int*)(aij->a + aij->nz); 2181 aij->i = aij->j + aij->nz; 2182 aij->singlemalloc = PETSC_TRUE; 2183 aij->freedata = PETSC_TRUE; 2184 2185 ierr = PetscMemcpy(aij->a,mxGetPr(mmat),aij->nz*sizeof(Scalar));CHKERRQ(ierr); 2186 /* Matlab stores by column, not row so we pass in the transpose of the matrix */ 2187 ierr = PetscMemcpy(aij->j,mxGetIr(mmat),aij->nz*sizeof(int));CHKERRQ(ierr); 2188 ierr = PetscMemcpy(aij->i,mxGetJc(mmat),(mat->m+1)*sizeof(int));CHKERRQ(ierr); 2189 2190 for (ii=0; ii<mat->m; ii++) { 2191 aij->ilen[ii] = aij->imax[ii] = aij->i[ii+1] - aij->i[ii]; 2192 } 2193 2194 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2195 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2196 2197 PetscFunctionReturn(0); 2198 } 2199 EXTERN_C_END 2200 #endif 2201 2202 /* --------------------------------------------------------------------------------*/ 2203 2204 #undef __FUNC__ 2205 #define __FUNC__ "MatCreateSeqAIJ" 2206 /*@C 2207 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 2208 (the default parallel PETSc format). For good matrix assembly performance 2209 the user should preallocate the matrix storage by setting the parameter nz 2210 (or the array nnz). By setting these parameters accurately, performance 2211 during matrix assembly can be increased by more than a factor of 50. 2212 2213 Collective on MPI_Comm 2214 2215 Input Parameters: 2216 + comm - MPI communicator, set to PETSC_COMM_SELF 2217 . m - number of rows 2218 . n - number of columns 2219 . nz - number of nonzeros per row (same for all rows) 2220 - nnz - array containing the number of nonzeros in the various rows 2221 (possibly different for each row) or PETSC_NULL 2222 2223 Output Parameter: 2224 . A - the matrix 2225 2226 Notes: 2227 The AIJ format (also called the Yale sparse matrix format or 2228 compressed row storage), is fully compatible with standard Fortran 77 2229 storage. That is, the stored row and column indices can begin at 2230 either one (as in Fortran) or zero. See the users' manual for details. 2231 2232 Specify the preallocated storage with either nz or nnz (not both). 2233 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2234 allocation. For large problems you MUST preallocate memory or you 2235 will get TERRIBLE performance, see the users' manual chapter on matrices. 2236 2237 By default, this format uses inodes (identical nodes) when possible, to 2238 improve numerical efficiency of matrix-vector products and solves. We 2239 search for consecutive rows with the same nonzero structure, thereby 2240 reusing matrix information to achieve increased efficiency. 2241 2242 Options Database Keys: 2243 + -mat_aij_no_inode - Do not use inodes 2244 . -mat_aij_inode_limit <limit> - Sets inode limit (max limit=5) 2245 - -mat_aij_oneindex - Internally use indexing starting at 1 2246 rather than 0. Note that when calling MatSetValues(), 2247 the user still MUST index entries starting at 0! 2248 2249 Level: intermediate 2250 2251 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2252 2253 @*/ 2254 int MatCreateSeqAIJ(MPI_Comm comm,int m,int n,int nz,int *nnz,Mat *A) 2255 { 2256 int ierr; 2257 2258 PetscFunctionBegin; 2259 ierr = MatCreate(comm,m,n,m,n,A);CHKERRQ(ierr); 2260 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 2261 ierr = MatSeqAIJSetPreallocation(*A,nz,nnz);CHKERRQ(ierr); 2262 PetscFunctionReturn(0); 2263 } 2264 2265 #undef __FUNC__ 2266 #define __FUNC__ "MatSeqAIJSetPreallocation" 2267 /*@C 2268 MatSeqAIJSetPreallocation - For good matrix assembly performance 2269 the user should preallocate the matrix storage by setting the parameter nz 2270 (or the array nnz). By setting these parameters accurately, performance 2271 during matrix assembly can be increased by more than a factor of 50. 2272 2273 Collective on MPI_Comm 2274 2275 Input Parameters: 2276 + comm - MPI communicator, set to PETSC_COMM_SELF 2277 . m - number of rows 2278 . n - number of columns 2279 . nz - number of nonzeros per row (same for all rows) 2280 - nnz - array containing the number of nonzeros in the various rows 2281 (possibly different for each row) or PETSC_NULL 2282 2283 Output Parameter: 2284 . A - the matrix 2285 2286 Notes: 2287 The AIJ format (also called the Yale sparse matrix format or 2288 compressed row storage), is fully compatible with standard Fortran 77 2289 storage. That is, the stored row and column indices can begin at 2290 either one (as in Fortran) or zero. See the users' manual for details. 2291 2292 Specify the preallocated storage with either nz or nnz (not both). 2293 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2294 allocation. For large problems you MUST preallocate memory or you 2295 will get TERRIBLE performance, see the users' manual chapter on matrices. 2296 2297 By default, this format uses inodes (identical nodes) when possible, to 2298 improve numerical efficiency of matrix-vector products and solves. We 2299 search for consecutive rows with the same nonzero structure, thereby 2300 reusing matrix information to achieve increased efficiency. 2301 2302 Options Database Keys: 2303 + -mat_aij_no_inode - Do not use inodes 2304 . -mat_aij_inode_limit <limit> - Sets inode limit (max limit=5) 2305 - -mat_aij_oneindex - Internally use indexing starting at 1 2306 rather than 0. Note that when calling MatSetValues(), 2307 the user still MUST index entries starting at 0! 2308 2309 Level: intermediate 2310 2311 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2312 2313 @*/ 2314 int MatSeqAIJSetPreallocation(Mat B,int nz,int *nnz) 2315 { 2316 Mat_SeqAIJ *b; 2317 int i,len,ierr; 2318 PetscTruth flg2; 2319 2320 PetscFunctionBegin; 2321 ierr = PetscTypeCompare((PetscObject)B,MATSEQAIJ,&flg2);CHKERRQ(ierr); 2322 if (!flg2) PetscFunctionReturn(0); 2323 2324 if (nz < -2) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than -2: value %d",nz); 2325 if (nnz) { 2326 for (i=0; i<B->m; i++) { 2327 if (nnz[i] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]); 2328 } 2329 } 2330 2331 B->preallocated = PETSC_TRUE; 2332 b = (Mat_SeqAIJ*)B->data; 2333 2334 ierr = PetscMalloc((B->m+1)*sizeof(int),&b->imax);CHKERRQ(ierr); 2335 if (!nnz) { 2336 if (nz == PETSC_DEFAULT) nz = 10; 2337 else if (nz <= 0) nz = 1; 2338 for (i=0; i<B->m; i++) b->imax[i] = nz; 2339 nz = nz*B->m; 2340 } else { 2341 nz = 0; 2342 for (i=0; i<B->m; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 2343 } 2344 2345 /* allocate the matrix space */ 2346 len = nz*(sizeof(int) + sizeof(Scalar)) + (B->m+1)*sizeof(int); 2347 ierr = PetscMalloc(len,&b->a);CHKERRQ(ierr); 2348 b->j = (int*)(b->a + nz); 2349 ierr = PetscMemzero(b->j,nz*sizeof(int));CHKERRQ(ierr); 2350 b->i = b->j + nz; 2351 b->singlemalloc = PETSC_TRUE; 2352 b->freedata = PETSC_TRUE; 2353 2354 b->i[0] = -b->indexshift; 2355 for (i=1; i<B->m+1; i++) { 2356 b->i[i] = b->i[i-1] + b->imax[i-1]; 2357 } 2358 2359 /* b->ilen will count nonzeros in each row so far. */ 2360 ierr = PetscMalloc((B->m+1)*sizeof(int),&b->ilen);CHKERRQ(ierr); 2361 PetscLogObjectMemory(B,len+2*(B->m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2362 for (i=0; i<B->m; i++) { b->ilen[i] = 0;} 2363 2364 b->nz = 0; 2365 b->maxnz = nz; 2366 B->info.nz_unneeded = (double)b->maxnz; 2367 PetscFunctionReturn(0); 2368 } 2369 2370 EXTERN_C_BEGIN 2371 #undef __FUNC__ 2372 #define __FUNC__ "MatCreate_SeqAIJ" 2373 int MatCreate_SeqAIJ(Mat B) 2374 { 2375 Mat_SeqAIJ *b; 2376 int ierr,size; 2377 PetscTruth flg; 2378 2379 PetscFunctionBegin; 2380 ierr = MPI_Comm_size(B->comm,&size);CHKERRQ(ierr); 2381 if (size > 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 2382 2383 B->m = B->M = PetscMax(B->m,B->M); 2384 B->n = B->N = PetscMax(B->n,B->N); 2385 2386 ierr = PetscNew(Mat_SeqAIJ,&b);CHKERRQ(ierr); 2387 B->data = (void*)b; 2388 ierr = PetscMemzero(b,sizeof(Mat_SeqAIJ));CHKERRQ(ierr); 2389 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2390 B->factor = 0; 2391 B->lupivotthreshold = 1.0; 2392 B->mapping = 0; 2393 ierr = PetscOptionsGetDouble(PETSC_NULL,"-mat_lu_pivotthreshold",&B->lupivotthreshold,PETSC_NULL);CHKERRQ(ierr); 2394 ierr = PetscOptionsHasName(PETSC_NULL,"-pc_ilu_preserve_row_sums",&b->ilu_preserve_row_sums);CHKERRQ(ierr); 2395 b->row = 0; 2396 b->col = 0; 2397 b->icol = 0; 2398 b->indexshift = 0; 2399 b->reallocs = 0; 2400 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_aij_oneindex",&flg);CHKERRQ(ierr); 2401 if (flg) b->indexshift = -1; 2402 2403 ierr = MapCreateMPI(B->comm,B->m,B->m,&B->rmap);CHKERRQ(ierr); 2404 ierr = MapCreateMPI(B->comm,B->n,B->n,&B->cmap);CHKERRQ(ierr); 2405 2406 b->sorted = PETSC_FALSE; 2407 b->ignorezeroentries = PETSC_FALSE; 2408 b->roworiented = PETSC_TRUE; 2409 b->nonew = 0; 2410 b->diag = 0; 2411 b->solve_work = 0; 2412 b->spptr = 0; 2413 b->inode.use = PETSC_TRUE; 2414 b->inode.node_count = 0; 2415 b->inode.size = 0; 2416 b->inode.limit = 5; 2417 b->inode.max_limit = 5; 2418 b->saved_values = 0; 2419 b->idiag = 0; 2420 b->ssor = 0; 2421 b->keepzeroedrows = PETSC_FALSE; 2422 2423 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 2424 2425 /* SuperLU is not currently supported through PETSc */ 2426 #if defined(PETSC_HAVE_SUPERLU) 2427 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_aij_superlu",&flg);CHKERRQ(ierr); 2428 if (flg) { ierr = MatUseSuperLU_SeqAIJ(B);CHKERRQ(ierr); } 2429 #endif 2430 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_aij_essl",&flg);CHKERRQ(ierr); 2431 if (flg) { ierr = MatUseEssl_SeqAIJ(B);CHKERRQ(ierr); } 2432 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_aij_lusol",&flg);CHKERRQ(ierr); 2433 if (flg) { ierr = MatUseLUSOL_SeqAIJ(B);CHKERRQ(ierr); } 2434 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_aij_matlab",&flg);CHKERRQ(ierr); 2435 if (flg) {ierr = MatUseMatlab_SeqAIJ(B);CHKERRQ(ierr);} 2436 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_aij_dxml",&flg);CHKERRQ(ierr); 2437 if (flg) { 2438 if (!b->indexshift) SETERRQ(PETSC_ERR_LIB,"need -mat_aij_oneindex with -mat_aij_dxml"); 2439 ierr = MatUseDXML_SeqAIJ(B);CHKERRQ(ierr); 2440 } 2441 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 2442 "MatSeqAIJSetColumnIndices_SeqAIJ", 2443 MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 2444 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 2445 "MatStoreValues_SeqAIJ", 2446 MatStoreValues_SeqAIJ);CHKERRQ(ierr); 2447 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 2448 "MatRetrieveValues_SeqAIJ", 2449 MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 2450 #if defined(PETSC_HAVE_MATLAB_ENGINE) && !defined(PETSC_USE_COMPLEX) 2451 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEnginePut_C","MatMatlabEnginePut_SeqAIJ",MatMatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 2452 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEngineGet_C","MatMatlabEngineGet_SeqAIJ",MatMatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 2453 #endif 2454 PetscFunctionReturn(0); 2455 } 2456 EXTERN_C_END 2457 2458 #undef __FUNC__ 2459 #define __FUNC__ "MatDuplicate_SeqAIJ" 2460 int MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 2461 { 2462 Mat C; 2463 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 2464 int i,len,m = A->m,shift = a->indexshift,ierr; 2465 2466 PetscFunctionBegin; 2467 *B = 0; 2468 ierr = MatCreate(A->comm,A->m,A->n,A->m,A->n,&C);CHKERRQ(ierr); 2469 ierr = MatSetType(C,MATSEQAIJ);CHKERRQ(ierr); 2470 c = (Mat_SeqAIJ*)C->data; 2471 2472 C->factor = A->factor; 2473 c->row = 0; 2474 c->col = 0; 2475 c->icol = 0; 2476 c->indexshift = shift; 2477 c->keepzeroedrows = a->keepzeroedrows; 2478 C->assembled = PETSC_TRUE; 2479 2480 C->M = A->m; 2481 C->N = A->n; 2482 2483 ierr = PetscMalloc((m+1)*sizeof(int),&c->imax);CHKERRQ(ierr); 2484 ierr = PetscMalloc((m+1)*sizeof(int),&c->ilen);CHKERRQ(ierr); 2485 for (i=0; i<m; i++) { 2486 c->imax[i] = a->imax[i]; 2487 c->ilen[i] = a->ilen[i]; 2488 } 2489 2490 /* allocate the matrix space */ 2491 c->singlemalloc = PETSC_TRUE; 2492 len = (m+1)*sizeof(int)+(a->i[m])*(sizeof(Scalar)+sizeof(int)); 2493 ierr = PetscMalloc(len,&c->a);CHKERRQ(ierr); 2494 c->j = (int*)(c->a + a->i[m] + shift); 2495 c->i = c->j + a->i[m] + shift; 2496 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(int));CHKERRQ(ierr); 2497 if (m > 0) { 2498 ierr = PetscMemcpy(c->j,a->j,(a->i[m]+shift)*sizeof(int));CHKERRQ(ierr); 2499 if (cpvalues == MAT_COPY_VALUES) { 2500 ierr = PetscMemcpy(c->a,a->a,(a->i[m]+shift)*sizeof(Scalar));CHKERRQ(ierr); 2501 } else { 2502 ierr = PetscMemzero(c->a,(a->i[m]+shift)*sizeof(Scalar));CHKERRQ(ierr); 2503 } 2504 } 2505 2506 PetscLogObjectMemory(C,len+2*(m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2507 c->sorted = a->sorted; 2508 c->roworiented = a->roworiented; 2509 c->nonew = a->nonew; 2510 c->ilu_preserve_row_sums = a->ilu_preserve_row_sums; 2511 c->saved_values = 0; 2512 c->idiag = 0; 2513 c->ssor = 0; 2514 c->ignorezeroentries = a->ignorezeroentries; 2515 c->freedata = PETSC_TRUE; 2516 2517 if (a->diag) { 2518 ierr = PetscMalloc((m+1)*sizeof(int),&c->diag);CHKERRQ(ierr); 2519 PetscLogObjectMemory(C,(m+1)*sizeof(int)); 2520 for (i=0; i<m; i++) { 2521 c->diag[i] = a->diag[i]; 2522 } 2523 } else c->diag = 0; 2524 c->inode.use = a->inode.use; 2525 c->inode.limit = a->inode.limit; 2526 c->inode.max_limit = a->inode.max_limit; 2527 if (a->inode.size){ 2528 ierr = PetscMalloc((m+1)*sizeof(int),&c->inode.size);CHKERRQ(ierr); 2529 c->inode.node_count = a->inode.node_count; 2530 ierr = PetscMemcpy(c->inode.size,a->inode.size,(m+1)*sizeof(int));CHKERRQ(ierr); 2531 } else { 2532 c->inode.size = 0; 2533 c->inode.node_count = 0; 2534 } 2535 c->nz = a->nz; 2536 c->maxnz = a->maxnz; 2537 c->solve_work = 0; 2538 c->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 2539 C->preallocated = PETSC_TRUE; 2540 2541 *B = C; 2542 ierr = PetscFListDuplicate(A->qlist,&C->qlist);CHKERRQ(ierr); 2543 PetscFunctionReturn(0); 2544 } 2545 2546 EXTERN_C_BEGIN 2547 #undef __FUNC__ 2548 #define __FUNC__ "MatLoad_SeqAIJ" 2549 int MatLoad_SeqAIJ(PetscViewer viewer,MatType type,Mat *A) 2550 { 2551 Mat_SeqAIJ *a; 2552 Mat B; 2553 int i,nz,ierr,fd,header[4],size,*rowlengths = 0,M,N,shift; 2554 MPI_Comm comm; 2555 2556 PetscFunctionBegin; 2557 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 2558 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2559 if (size > 1) SETERRQ(PETSC_ERR_ARG_SIZ,"view must have one processor"); 2560 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 2561 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 2562 if (header[0] != MAT_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 2563 M = header[1]; N = header[2]; nz = header[3]; 2564 2565 if (nz < 0) { 2566 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 2567 } 2568 2569 /* read in row lengths */ 2570 ierr = PetscMalloc(M*sizeof(int),&rowlengths);CHKERRQ(ierr); 2571 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 2572 2573 /* create our matrix */ 2574 ierr = MatCreateSeqAIJ(comm,M,N,0,rowlengths,A);CHKERRQ(ierr); 2575 B = *A; 2576 a = (Mat_SeqAIJ*)B->data; 2577 shift = a->indexshift; 2578 2579 /* read in column indices and adjust for Fortran indexing*/ 2580 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 2581 if (shift) { 2582 for (i=0; i<nz; i++) { 2583 a->j[i] += 1; 2584 } 2585 } 2586 2587 /* read in nonzero values */ 2588 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 2589 2590 /* set matrix "i" values */ 2591 a->i[0] = -shift; 2592 for (i=1; i<= M; i++) { 2593 a->i[i] = a->i[i-1] + rowlengths[i-1]; 2594 a->ilen[i-1] = rowlengths[i-1]; 2595 } 2596 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 2597 2598 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2599 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2600 PetscFunctionReturn(0); 2601 } 2602 EXTERN_C_END 2603 2604 #undef __FUNC__ 2605 #define __FUNC__ /*<a name="MatEqual_SeqAIJ""></a>*/"MatEqual_SeqAIJ" 2606 int MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg) 2607 { 2608 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data,*b = (Mat_SeqAIJ *)B->data; 2609 int ierr; 2610 PetscTruth flag; 2611 2612 PetscFunctionBegin; 2613 ierr = PetscTypeCompare((PetscObject)B,MATSEQAIJ,&flag);CHKERRQ(ierr); 2614 if (!flag) SETERRQ(PETSC_ERR_ARG_INCOMP,"Matrices must be same type"); 2615 2616 /* If the matrix dimensions are not equal,or no of nonzeros or shift */ 2617 if ((A->m != B->m) || (A->n != B->n) ||(a->nz != b->nz)|| (a->indexshift != b->indexshift)) { 2618 *flg = PETSC_FALSE; 2619 PetscFunctionReturn(0); 2620 } 2621 2622 /* if the a->i are the same */ 2623 ierr = PetscMemcmp(a->i,b->i,(A->m+1)*sizeof(int),flg);CHKERRQ(ierr); 2624 if (*flg == PETSC_FALSE) PetscFunctionReturn(0); 2625 2626 /* if a->j are the same */ 2627 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(int),flg);CHKERRQ(ierr); 2628 if (*flg == PETSC_FALSE) PetscFunctionReturn(0); 2629 2630 /* if a->a are the same */ 2631 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(Scalar),flg);CHKERRQ(ierr); 2632 2633 PetscFunctionReturn(0); 2634 2635 } 2636 2637 #undef __FUNC__ 2638 #define __FUNC__ "MatCreateSeqAIJWithArrays" 2639 /*@C 2640 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 2641 provided by the user. 2642 2643 Coolective on MPI_Comm 2644 2645 Input Parameters: 2646 + comm - must be an MPI communicator of size 1 2647 . m - number of rows 2648 . n - number of columns 2649 . i - row indices 2650 . j - column indices 2651 - a - matrix values 2652 2653 Output Parameter: 2654 . mat - the matrix 2655 2656 Level: intermediate 2657 2658 Notes: 2659 The i, j, and a arrays are not copied by this routine, the user must free these arrays 2660 once the matrix is destroyed 2661 2662 You cannot set new nonzero locations into this matrix, that will generate an error. 2663 2664 The i and j indices can be either 0- or 1 based 2665 2666 .seealso: MatCreate(), MatCreateMPIAIJ(), MatCreateSeqAIJ() 2667 2668 @*/ 2669 int MatCreateSeqAIJWithArrays(MPI_Comm comm,int m,int n,int* i,int*j,Scalar *a,Mat *mat) 2670 { 2671 int ierr,ii; 2672 Mat_SeqAIJ *aij; 2673 2674 PetscFunctionBegin; 2675 ierr = MatCreateSeqAIJ(comm,m,n,0,0,mat);CHKERRQ(ierr); 2676 aij = (Mat_SeqAIJ*)(*mat)->data; 2677 ierr = PetscFree(aij->a);CHKERRQ(ierr); 2678 2679 if (i[0] == 1) { 2680 aij->indexshift = -1; 2681 } else if (i[0]) { 2682 SETERRQ(1,"i (row indices) do not start with 0 or 1"); 2683 } 2684 aij->i = i; 2685 aij->j = j; 2686 aij->a = a; 2687 aij->singlemalloc = PETSC_FALSE; 2688 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 2689 aij->freedata = PETSC_FALSE; 2690 2691 for (ii=0; ii<m; ii++) { 2692 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 2693 #if defined(PETSC_USE_BOPT_g) 2694 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]); 2695 #endif 2696 } 2697 #if defined(PETSC_USE_BOPT_g) 2698 for (ii=0; ii<aij->i[m]; ii++) { 2699 if (j[ii] < -aij->indexshift) SETERRQ2(1,"Negative column index at location = %d index = %d",ii,j[ii]); 2700 if (j[ii] > n - 1 -aij->indexshift) SETERRQ2(1,"Column index to large at location = %d index = %d",ii,j[ii]); 2701 } 2702 #endif 2703 2704 /* changes indices to start at 0 */ 2705 if (i[0]) { 2706 aij->indexshift = 0; 2707 for (ii=0; ii<m; ii++) { 2708 i[ii]--; 2709 } 2710 for (ii=0; ii<i[m]; ii++) { 2711 j[ii]--; 2712 } 2713 } 2714 2715 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2716 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2717 PetscFunctionReturn(0); 2718 } 2719 2720 2721 2722