1 /*$Id: aij.c,v 1.384 2001/08/29 16:12:22 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" /*I "petscmat.h" I*/ 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 __FUNCT__ 18 #define __FUNCT__ "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] - 1; 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]; 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 __FUNCT__ 51 #define __FUNCT__ "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 __FUNCT__ 67 #define __FUNCT__ "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 __FUNCT__ 108 #define __FUNCT__ "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 __FUNCT__ 125 #define __FUNCT__ "MatSetValues_SeqAIJ" 126 int MatSetValues_SeqAIJ(Mat A,int m,int *im,int n,int *in,PetscScalar *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 PetscScalar *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 PetscScalar *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(PetscScalar))+(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(PetscScalar));CHKERRQ(ierr); 195 ierr = PetscMemcpy(new_a+ai[row]+shift+nrow+CHUNKSIZE,aa+ai[row]+shift+nrow,len*sizeof(PetscScalar));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(PetscScalar))); 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 __FUNCT__ 228 #define __FUNCT__ "MatGetValues_SeqAIJ" 229 int MatGetValues_SeqAIJ(Mat A,int m,int *im,int n,int *in,PetscScalar *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 PetscScalar *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 __FUNCT__ 269 #define __FUNCT__ "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 __FUNCT__ 305 #define __FUNCT__ "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)A,&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 PetscScalar 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 __FUNCT__ 466 #define __FUNCT__ "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 __FUNCT__ 548 #define __FUNCT__ "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 __FUNCT__ 571 #define __FUNCT__ "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,PetscTruth); 601 #undef __FUNCT__ 602 #define __FUNCT__ "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 PetscScalar *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,(PetscTruth)(!fshift));CHKERRQ(ierr); 654 PetscFunctionReturn(0); 655 } 656 657 #undef __FUNCT__ 658 #define __FUNCT__ "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(PetscScalar));CHKERRQ(ierr); 666 PetscFunctionReturn(0); 667 } 668 669 #undef __FUNCT__ 670 #define __FUNCT__ "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 if (a->coloring) {ierr = ISColoringDestroy(a->coloring);CHKERRQ(ierr);} 702 ierr = PetscFree(a);CHKERRQ(ierr); 703 PetscFunctionReturn(0); 704 } 705 706 #undef __FUNCT__ 707 #define __FUNCT__ "MatCompress_SeqAIJ" 708 int MatCompress_SeqAIJ(Mat A) 709 { 710 PetscFunctionBegin; 711 PetscFunctionReturn(0); 712 } 713 714 #undef __FUNCT__ 715 #define __FUNCT__ "MatSetOption_SeqAIJ" 716 int MatSetOption_SeqAIJ(Mat A,MatOption op) 717 { 718 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 719 720 PetscFunctionBegin; 721 switch (op) { 722 case MAT_ROW_ORIENTED: 723 a->roworiented = PETSC_TRUE; 724 break; 725 case MAT_KEEP_ZEROED_ROWS: 726 a->keepzeroedrows = PETSC_TRUE; 727 break; 728 case MAT_COLUMN_ORIENTED: 729 a->roworiented = PETSC_FALSE; 730 break; 731 case MAT_COLUMNS_SORTED: 732 a->sorted = PETSC_TRUE; 733 break; 734 case MAT_COLUMNS_UNSORTED: 735 a->sorted = PETSC_FALSE; 736 break; 737 case MAT_NO_NEW_NONZERO_LOCATIONS: 738 a->nonew = 1; 739 break; 740 case MAT_NEW_NONZERO_LOCATION_ERR: 741 a->nonew = -1; 742 break; 743 case MAT_NEW_NONZERO_ALLOCATION_ERR: 744 a->nonew = -2; 745 break; 746 case MAT_YES_NEW_NONZERO_LOCATIONS: 747 a->nonew = 0; 748 break; 749 case MAT_IGNORE_ZERO_ENTRIES: 750 a->ignorezeroentries = PETSC_TRUE; 751 break; 752 case MAT_USE_INODES: 753 a->inode.use = PETSC_TRUE; 754 break; 755 case MAT_DO_NOT_USE_INODES: 756 a->inode.use = PETSC_FALSE; 757 break; 758 case MAT_ROWS_SORTED: 759 case MAT_ROWS_UNSORTED: 760 case MAT_YES_NEW_DIAGONALS: 761 case MAT_IGNORE_OFF_PROC_ENTRIES: 762 case MAT_USE_HASH_TABLE: 763 case MAT_USE_SINGLE_PRECISION_SOLVES: 764 PetscLogInfo(A,"MatSetOption_SeqAIJ:Option ignored\n"); 765 break; 766 case MAT_NO_NEW_DIAGONALS: 767 SETERRQ(PETSC_ERR_SUP,"MAT_NO_NEW_DIAGONALS"); 768 break; 769 case MAT_INODE_LIMIT_1: 770 a->inode.limit = 1; 771 break; 772 case MAT_INODE_LIMIT_2: 773 a->inode.limit = 2; 774 break; 775 case MAT_INODE_LIMIT_3: 776 a->inode.limit = 3; 777 break; 778 case MAT_INODE_LIMIT_4: 779 a->inode.limit = 4; 780 break; 781 case MAT_INODE_LIMIT_5: 782 a->inode.limit = 5; 783 break; 784 default: 785 SETERRQ(PETSC_ERR_SUP,"unknown option"); 786 break; 787 } 788 PetscFunctionReturn(0); 789 } 790 791 #undef __FUNCT__ 792 #define __FUNCT__ "MatGetDiagonal_SeqAIJ" 793 int MatGetDiagonal_SeqAIJ(Mat A,Vec v) 794 { 795 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 796 int i,j,n,shift = a->indexshift,ierr; 797 PetscScalar *x,zero = 0.0; 798 799 PetscFunctionBegin; 800 ierr = VecSet(&zero,v);CHKERRQ(ierr); 801 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 802 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 803 if (n != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 804 for (i=0; i<A->m; i++) { 805 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 806 if (a->j[j]+shift == i) { 807 x[i] = a->a[j]; 808 break; 809 } 810 } 811 } 812 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 813 PetscFunctionReturn(0); 814 } 815 816 817 #undef __FUNCT__ 818 #define __FUNCT__ "MatMultTransposeAdd_SeqAIJ" 819 int MatMultTransposeAdd_SeqAIJ(Mat A,Vec xx,Vec zz,Vec yy) 820 { 821 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 822 PetscScalar *x,*y; 823 int ierr,m = A->m,shift = a->indexshift; 824 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 825 PetscScalar *v,alpha; 826 int n,i,*idx; 827 #endif 828 829 PetscFunctionBegin; 830 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 831 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 832 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 833 y = y + shift; /* shift for Fortran start by 1 indexing */ 834 835 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 836 fortranmulttransposeaddaij_(&m,x,a->i,a->j+shift,a->a+shift,y); 837 #else 838 for (i=0; i<m; i++) { 839 idx = a->j + a->i[i] + shift; 840 v = a->a + a->i[i] + shift; 841 n = a->i[i+1] - a->i[i]; 842 alpha = x[i]; 843 while (n-->0) {y[*idx++] += alpha * *v++;} 844 } 845 #endif 846 PetscLogFlops(2*a->nz); 847 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 848 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 849 PetscFunctionReturn(0); 850 } 851 852 #undef __FUNCT__ 853 #define __FUNCT__ "MatMultTranspose_SeqAIJ" 854 int MatMultTranspose_SeqAIJ(Mat A,Vec xx,Vec yy) 855 { 856 PetscScalar zero = 0.0; 857 int ierr; 858 859 PetscFunctionBegin; 860 ierr = VecSet(&zero,yy);CHKERRQ(ierr); 861 ierr = MatMultTransposeAdd_SeqAIJ(A,xx,yy,yy);CHKERRQ(ierr); 862 PetscFunctionReturn(0); 863 } 864 865 866 #undef __FUNCT__ 867 #define __FUNCT__ "MatMult_SeqAIJ" 868 int MatMult_SeqAIJ(Mat A,Vec xx,Vec yy) 869 { 870 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 871 PetscScalar *x,*y,*v,sum; 872 int ierr,m = A->m,*idx,shift = a->indexshift,*ii; 873 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 874 int n,i,jrow,j; 875 #endif 876 877 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 878 #pragma disjoint(*x,*y,*v) 879 #endif 880 881 PetscFunctionBegin; 882 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 883 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 884 x = x + shift; /* shift for Fortran start by 1 indexing */ 885 idx = a->j; 886 v = a->a; 887 ii = a->i; 888 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 889 fortranmultaij_(&m,x,ii,idx+shift,v+shift,y); 890 #else 891 v += shift; /* shift for Fortran start by 1 indexing */ 892 idx += shift; 893 for (i=0; i<m; i++) { 894 jrow = ii[i]; 895 n = ii[i+1] - jrow; 896 sum = 0.0; 897 for (j=0; j<n; j++) { 898 sum += v[jrow]*x[idx[jrow]]; jrow++; 899 } 900 y[i] = sum; 901 } 902 #endif 903 PetscLogFlops(2*a->nz - m); 904 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 905 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 906 PetscFunctionReturn(0); 907 } 908 909 #undef __FUNCT__ 910 #define __FUNCT__ "MatMultAdd_SeqAIJ" 911 int MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 912 { 913 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 914 PetscScalar *x,*y,*z,*v,sum; 915 int ierr,m = A->m,*idx,shift = a->indexshift,*ii; 916 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 917 int n,i,jrow,j; 918 #endif 919 920 PetscFunctionBegin; 921 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 922 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 923 if (zz != yy) { 924 ierr = VecGetArray(zz,&z);CHKERRQ(ierr); 925 } else { 926 z = y; 927 } 928 x = x + shift; /* shift for Fortran start by 1 indexing */ 929 idx = a->j; 930 v = a->a; 931 ii = a->i; 932 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 933 fortranmultaddaij_(&m,x,ii,idx+shift,v+shift,y,z); 934 #else 935 v += shift; /* shift for Fortran start by 1 indexing */ 936 idx += shift; 937 for (i=0; i<m; i++) { 938 jrow = ii[i]; 939 n = ii[i+1] - jrow; 940 sum = y[i]; 941 for (j=0; j<n; j++) { 942 sum += v[jrow]*x[idx[jrow]]; jrow++; 943 } 944 z[i] = sum; 945 } 946 #endif 947 PetscLogFlops(2*a->nz); 948 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 949 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 950 if (zz != yy) { 951 ierr = VecRestoreArray(zz,&z);CHKERRQ(ierr); 952 } 953 PetscFunctionReturn(0); 954 } 955 956 /* 957 Adds diagonal pointers to sparse matrix structure. 958 */ 959 #undef __FUNCT__ 960 #define __FUNCT__ "MatMarkDiagonal_SeqAIJ" 961 int MatMarkDiagonal_SeqAIJ(Mat A) 962 { 963 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 964 int i,j,*diag,m = A->m,shift = a->indexshift,ierr; 965 966 PetscFunctionBegin; 967 if (a->diag) PetscFunctionReturn(0); 968 969 ierr = PetscMalloc((m+1)*sizeof(int),&diag);CHKERRQ(ierr); 970 PetscLogObjectMemory(A,(m+1)*sizeof(int)); 971 for (i=0; i<A->m; i++) { 972 diag[i] = a->i[i+1]; 973 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 974 if (a->j[j]+shift == i) { 975 diag[i] = j - shift; 976 break; 977 } 978 } 979 } 980 a->diag = diag; 981 PetscFunctionReturn(0); 982 } 983 984 /* 985 Checks for missing diagonals 986 */ 987 #undef __FUNCT__ 988 #define __FUNCT__ "MatMissingDiagonal_SeqAIJ" 989 int MatMissingDiagonal_SeqAIJ(Mat A) 990 { 991 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 992 int *diag,*jj = a->j,i,shift = a->indexshift,ierr; 993 994 PetscFunctionBegin; 995 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 996 diag = a->diag; 997 for (i=0; i<A->m; i++) { 998 if (jj[diag[i]+shift] != i-shift) { 999 SETERRQ1(1,"Matrix is missing diagonal number %d",i); 1000 } 1001 } 1002 PetscFunctionReturn(0); 1003 } 1004 1005 #undef __FUNCT__ 1006 #define __FUNCT__ "MatRelax_SeqAIJ" 1007 int MatRelax_SeqAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,int its,Vec xx) 1008 { 1009 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1010 PetscScalar *x,*b,*bs, d,*xs,sum,*v = a->a,*t=0,scale,*ts,*xb,*idiag=0; 1011 int ierr,*idx,*diag,n = A->n,m = A->m,i,shift = a->indexshift; 1012 1013 PetscFunctionBegin; 1014 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1015 if (xx != bb) { 1016 ierr = VecGetArray(bb,&b);CHKERRQ(ierr); 1017 } else { 1018 b = x; 1019 } 1020 1021 if (!a->diag) {ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr);} 1022 diag = a->diag; 1023 xs = x + shift; /* shifted by one for index start of a or a->j*/ 1024 if (flag == SOR_APPLY_UPPER) { 1025 /* apply (U + D/omega) to the vector */ 1026 bs = b + shift; 1027 for (i=0; i<m; i++) { 1028 d = fshift + a->a[diag[i] + shift]; 1029 n = a->i[i+1] - diag[i] - 1; 1030 PetscLogFlops(2*n-1); 1031 idx = a->j + diag[i] + (!shift); 1032 v = a->a + diag[i] + (!shift); 1033 sum = b[i]*d/omega; 1034 SPARSEDENSEDOT(sum,bs,v,idx,n); 1035 x[i] = sum; 1036 } 1037 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1038 if (bb != xx) {ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);} 1039 PetscFunctionReturn(0); 1040 } 1041 1042 /* setup workspace for Eisenstat */ 1043 if (flag & SOR_EISENSTAT) { 1044 if (!a->idiag) { 1045 ierr = PetscMalloc(2*m*sizeof(PetscScalar),&a->idiag);CHKERRQ(ierr); 1046 a->ssor = a->idiag + m; 1047 v = a->a; 1048 for (i=0; i<m; i++) { a->idiag[i] = 1.0/v[diag[i]];} 1049 } 1050 t = a->ssor; 1051 idiag = a->idiag; 1052 } 1053 /* Let A = L + U + D; where L is lower trianglar, 1054 U is upper triangular, E is diagonal; This routine applies 1055 1056 (L + E)^{-1} A (U + E)^{-1} 1057 1058 to a vector efficiently using Eisenstat's trick. This is for 1059 the case of SSOR preconditioner, so E is D/omega where omega 1060 is the relaxation factor. 1061 */ 1062 1063 if (flag == SOR_APPLY_LOWER) { 1064 SETERRQ(PETSC_ERR_SUP,"SOR_APPLY_LOWER is not implemented"); 1065 } else if ((flag & SOR_EISENSTAT) && omega == 1.0 && shift == 0 && fshift == 0.0) { 1066 /* special case for omega = 1.0 saves flops and some integer ops */ 1067 PetscScalar *v2; 1068 1069 v2 = a->a; 1070 /* x = (E + U)^{-1} b */ 1071 for (i=m-1; i>=0; i--) { 1072 n = a->i[i+1] - diag[i] - 1; 1073 idx = a->j + diag[i] + 1; 1074 v = a->a + diag[i] + 1; 1075 sum = b[i]; 1076 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1077 x[i] = sum*idiag[i]; 1078 1079 /* t = b - (2*E - D)x */ 1080 t[i] = b[i] - (v2[diag[i]])*x[i]; 1081 } 1082 1083 /* t = (E + L)^{-1}t */ 1084 diag = a->diag; 1085 for (i=0; i<m; i++) { 1086 n = diag[i] - a->i[i]; 1087 idx = a->j + a->i[i]; 1088 v = a->a + a->i[i]; 1089 sum = t[i]; 1090 SPARSEDENSEMDOT(sum,t,v,idx,n); 1091 t[i] = sum*idiag[i]; 1092 1093 /* x = x + t */ 1094 x[i] += t[i]; 1095 } 1096 1097 PetscLogFlops(3*m-1 + 2*a->nz); 1098 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1099 if (bb != xx) {ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);} 1100 PetscFunctionReturn(0); 1101 } else if (flag & SOR_EISENSTAT) { 1102 /* Let A = L + U + D; where L is lower trianglar, 1103 U is upper triangular, E is diagonal; This routine applies 1104 1105 (L + E)^{-1} A (U + E)^{-1} 1106 1107 to a vector efficiently using Eisenstat's trick. This is for 1108 the case of SSOR preconditioner, so E is D/omega where omega 1109 is the relaxation factor. 1110 */ 1111 scale = (2.0/omega) - 1.0; 1112 1113 /* x = (E + U)^{-1} b */ 1114 for (i=m-1; i>=0; i--) { 1115 d = fshift + a->a[diag[i] + shift]; 1116 n = a->i[i+1] - diag[i] - 1; 1117 idx = a->j + diag[i] + (!shift); 1118 v = a->a + diag[i] + (!shift); 1119 sum = b[i]; 1120 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1121 x[i] = omega*(sum/d); 1122 } 1123 1124 /* t = b - (2*E - D)x */ 1125 v = a->a; 1126 for (i=0; i<m; i++) { t[i] = b[i] - scale*(v[*diag++ + shift])*x[i]; } 1127 1128 /* t = (E + L)^{-1}t */ 1129 ts = t + shift; /* shifted by one for index start of a or a->j*/ 1130 diag = a->diag; 1131 for (i=0; i<m; i++) { 1132 d = fshift + a->a[diag[i]+shift]; 1133 n = diag[i] - a->i[i]; 1134 idx = a->j + a->i[i] + shift; 1135 v = a->a + a->i[i] + shift; 1136 sum = t[i]; 1137 SPARSEDENSEMDOT(sum,ts,v,idx,n); 1138 t[i] = omega*(sum/d); 1139 /* x = x + t */ 1140 x[i] += t[i]; 1141 } 1142 1143 PetscLogFlops(6*m-1 + 2*a->nz); 1144 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1145 if (bb != xx) {ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);} 1146 PetscFunctionReturn(0); 1147 } 1148 if (flag & SOR_ZERO_INITIAL_GUESS) { 1149 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1150 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1151 fortranrelaxaijforwardzero_(&m,&omega,x,a->i,a->j,diag,a->a,b); 1152 #else 1153 for (i=0; i<m; i++) { 1154 d = fshift + a->a[diag[i]+shift]; 1155 n = diag[i] - a->i[i]; 1156 PetscLogFlops(2*n-1); 1157 idx = a->j + a->i[i] + shift; 1158 v = a->a + a->i[i] + shift; 1159 sum = b[i]; 1160 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1161 x[i] = omega*(sum/d); 1162 } 1163 #endif 1164 xb = x; 1165 } else xb = b; 1166 if ((flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) && 1167 (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP)) { 1168 for (i=0; i<m; i++) { 1169 x[i] *= a->a[diag[i]+shift]; 1170 } 1171 PetscLogFlops(m); 1172 } 1173 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1174 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1175 fortranrelaxaijbackwardzero_(&m,&omega,x,a->i,a->j,diag,a->a,xb); 1176 #else 1177 for (i=m-1; i>=0; i--) { 1178 d = fshift + a->a[diag[i] + shift]; 1179 n = a->i[i+1] - diag[i] - 1; 1180 PetscLogFlops(2*n-1); 1181 idx = a->j + diag[i] + (!shift); 1182 v = a->a + diag[i] + (!shift); 1183 sum = xb[i]; 1184 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1185 x[i] = omega*(sum/d); 1186 } 1187 #endif 1188 } 1189 its--; 1190 } 1191 while (its--) { 1192 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1193 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1194 fortranrelaxaijforward_(&m,&omega,x,a->i,a->j,diag,a->a,b); 1195 #else 1196 for (i=0; i<m; i++) { 1197 d = fshift + a->a[diag[i]+shift]; 1198 n = a->i[i+1] - a->i[i]; 1199 PetscLogFlops(2*n-1); 1200 idx = a->j + a->i[i] + shift; 1201 v = a->a + a->i[i] + shift; 1202 sum = b[i]; 1203 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1204 x[i] = (1. - omega)*x[i] + omega*(sum + a->a[diag[i]+shift]*x[i])/d; 1205 } 1206 #endif 1207 } 1208 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1209 #if defined(PETSC_USE_FORTRAN_KERNEL_RELAXAIJ) 1210 fortranrelaxaijbackward_(&m,&omega,x,a->i,a->j,diag,a->a,b); 1211 #else 1212 for (i=m-1; i>=0; i--) { 1213 d = fshift + a->a[diag[i] + shift]; 1214 n = a->i[i+1] - a->i[i]; 1215 PetscLogFlops(2*n-1); 1216 idx = a->j + a->i[i] + shift; 1217 v = a->a + a->i[i] + shift; 1218 sum = b[i]; 1219 SPARSEDENSEMDOT(sum,xs,v,idx,n); 1220 x[i] = (1. - omega)*x[i] + omega*(sum + a->a[diag[i]+shift]*x[i])/d; 1221 } 1222 #endif 1223 } 1224 } 1225 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1226 if (bb != xx) {ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr);} 1227 PetscFunctionReturn(0); 1228 } 1229 1230 #undef __FUNCT__ 1231 #define __FUNCT__ "MatGetInfo_SeqAIJ" 1232 int MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1233 { 1234 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1235 1236 PetscFunctionBegin; 1237 info->rows_global = (double)A->m; 1238 info->columns_global = (double)A->n; 1239 info->rows_local = (double)A->m; 1240 info->columns_local = (double)A->n; 1241 info->block_size = 1.0; 1242 info->nz_allocated = (double)a->maxnz; 1243 info->nz_used = (double)a->nz; 1244 info->nz_unneeded = (double)(a->maxnz - a->nz); 1245 info->assemblies = (double)A->num_ass; 1246 info->mallocs = (double)a->reallocs; 1247 info->memory = A->mem; 1248 if (A->factor) { 1249 info->fill_ratio_given = A->info.fill_ratio_given; 1250 info->fill_ratio_needed = A->info.fill_ratio_needed; 1251 info->factor_mallocs = A->info.factor_mallocs; 1252 } else { 1253 info->fill_ratio_given = 0; 1254 info->fill_ratio_needed = 0; 1255 info->factor_mallocs = 0; 1256 } 1257 PetscFunctionReturn(0); 1258 } 1259 1260 EXTERN int MatLUFactorSymbolic_SeqAIJ(Mat,IS,IS,MatLUInfo*,Mat*); 1261 EXTERN int MatLUFactorNumeric_SeqAIJ(Mat,Mat*); 1262 EXTERN int MatLUFactor_SeqAIJ(Mat,IS,IS,MatLUInfo*); 1263 EXTERN int MatSolve_SeqAIJ(Mat,Vec,Vec); 1264 EXTERN int MatSolveAdd_SeqAIJ(Mat,Vec,Vec,Vec); 1265 EXTERN int MatSolveTranspose_SeqAIJ(Mat,Vec,Vec); 1266 EXTERN int MatSolveTransposeAdd_SeqAIJ(Mat,Vec,Vec,Vec); 1267 1268 #undef __FUNCT__ 1269 #define __FUNCT__ "MatZeroRows_SeqAIJ" 1270 int MatZeroRows_SeqAIJ(Mat A,IS is,PetscScalar *diag) 1271 { 1272 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1273 int i,ierr,N,*rows,m = A->m - 1,shift = a->indexshift; 1274 1275 PetscFunctionBegin; 1276 ierr = ISGetLocalSize(is,&N);CHKERRQ(ierr); 1277 ierr = ISGetIndices(is,&rows);CHKERRQ(ierr); 1278 if (a->keepzeroedrows) { 1279 for (i=0; i<N; i++) { 1280 if (rows[i] < 0 || rows[i] > m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"row out of range"); 1281 ierr = PetscMemzero(&a->a[a->i[rows[i]]+shift],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 1282 } 1283 if (diag) { 1284 ierr = MatMissingDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1285 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1286 for (i=0; i<N; i++) { 1287 a->a[a->diag[rows[i]]] = *diag; 1288 } 1289 } 1290 } else { 1291 if (diag) { 1292 for (i=0; i<N; i++) { 1293 if (rows[i] < 0 || rows[i] > m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"row out of range"); 1294 if (a->ilen[rows[i]] > 0) { 1295 a->ilen[rows[i]] = 1; 1296 a->a[a->i[rows[i]]+shift] = *diag; 1297 a->j[a->i[rows[i]]+shift] = rows[i]+shift; 1298 } else { /* in case row was completely empty */ 1299 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],diag,INSERT_VALUES);CHKERRQ(ierr); 1300 } 1301 } 1302 } else { 1303 for (i=0; i<N; i++) { 1304 if (rows[i] < 0 || rows[i] > m) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"row out of range"); 1305 a->ilen[rows[i]] = 0; 1306 } 1307 } 1308 } 1309 ierr = ISRestoreIndices(is,&rows);CHKERRQ(ierr); 1310 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1311 PetscFunctionReturn(0); 1312 } 1313 1314 #undef __FUNCT__ 1315 #define __FUNCT__ "MatGetOwnershipRange_SeqAIJ" 1316 int MatGetOwnershipRange_SeqAIJ(Mat A,int *m,int *n) 1317 { 1318 PetscFunctionBegin; 1319 if (m) *m = 0; 1320 if (n) *n = A->m; 1321 PetscFunctionReturn(0); 1322 } 1323 1324 #undef __FUNCT__ 1325 #define __FUNCT__ "MatGetRow_SeqAIJ" 1326 int MatGetRow_SeqAIJ(Mat A,int row,int *nz,int **idx,PetscScalar **v) 1327 { 1328 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1329 int *itmp,i,shift = a->indexshift,ierr; 1330 1331 PetscFunctionBegin; 1332 if (row < 0 || row >= A->m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Row %d out of range",row); 1333 1334 *nz = a->i[row+1] - a->i[row]; 1335 if (v) *v = a->a + a->i[row] + shift; 1336 if (idx) { 1337 itmp = a->j + a->i[row] + shift; 1338 if (*nz && shift) { 1339 ierr = PetscMalloc((*nz)*sizeof(int),idx);CHKERRQ(ierr); 1340 for (i=0; i<(*nz); i++) {(*idx)[i] = itmp[i] + shift;} 1341 } else if (*nz) { 1342 *idx = itmp; 1343 } 1344 else *idx = 0; 1345 } 1346 PetscFunctionReturn(0); 1347 } 1348 1349 #undef __FUNCT__ 1350 #define __FUNCT__ "MatRestoreRow_SeqAIJ" 1351 int MatRestoreRow_SeqAIJ(Mat A,int row,int *nz,int **idx,PetscScalar **v) 1352 { 1353 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1354 int ierr; 1355 1356 PetscFunctionBegin; 1357 if (idx) {if (*idx && a->indexshift) {ierr = PetscFree(*idx);CHKERRQ(ierr);}} 1358 PetscFunctionReturn(0); 1359 } 1360 1361 #undef __FUNCT__ 1362 #define __FUNCT__ "MatNorm_SeqAIJ" 1363 int MatNorm_SeqAIJ(Mat A,NormType type,PetscReal *norm) 1364 { 1365 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1366 PetscScalar *v = a->a; 1367 PetscReal sum = 0.0; 1368 int i,j,shift = a->indexshift,ierr; 1369 1370 PetscFunctionBegin; 1371 if (type == NORM_FROBENIUS) { 1372 for (i=0; i<a->nz; i++) { 1373 #if defined(PETSC_USE_COMPLEX) 1374 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1375 #else 1376 sum += (*v)*(*v); v++; 1377 #endif 1378 } 1379 *norm = sqrt(sum); 1380 } else if (type == NORM_1) { 1381 PetscReal *tmp; 1382 int *jj = a->j; 1383 ierr = PetscMalloc((A->n+1)*sizeof(PetscReal),&tmp);CHKERRQ(ierr); 1384 ierr = PetscMemzero(tmp,A->n*sizeof(PetscReal));CHKERRQ(ierr); 1385 *norm = 0.0; 1386 for (j=0; j<a->nz; j++) { 1387 tmp[*jj++ + shift] += PetscAbsScalar(*v); v++; 1388 } 1389 for (j=0; j<A->n; j++) { 1390 if (tmp[j] > *norm) *norm = tmp[j]; 1391 } 1392 ierr = PetscFree(tmp);CHKERRQ(ierr); 1393 } else if (type == NORM_INFINITY) { 1394 *norm = 0.0; 1395 for (j=0; j<A->m; j++) { 1396 v = a->a + a->i[j] + shift; 1397 sum = 0.0; 1398 for (i=0; i<a->i[j+1]-a->i[j]; i++) { 1399 sum += PetscAbsScalar(*v); v++; 1400 } 1401 if (sum > *norm) *norm = sum; 1402 } 1403 } else { 1404 SETERRQ(PETSC_ERR_SUP,"No support for two norm"); 1405 } 1406 PetscFunctionReturn(0); 1407 } 1408 1409 #undef __FUNCT__ 1410 #define __FUNCT__ "MatTranspose_SeqAIJ" 1411 int MatTranspose_SeqAIJ(Mat A,Mat *B) 1412 { 1413 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1414 Mat C; 1415 int i,ierr,*aj = a->j,*ai = a->i,m = A->m,len,*col; 1416 int shift = a->indexshift; 1417 PetscScalar *array = a->a; 1418 1419 PetscFunctionBegin; 1420 if (!B && m != A->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Square matrix only for in-place"); 1421 ierr = PetscMalloc((1+A->n)*sizeof(int),&col);CHKERRQ(ierr); 1422 ierr = PetscMemzero(col,(1+A->n)*sizeof(int));CHKERRQ(ierr); 1423 if (shift) { 1424 for (i=0; i<ai[m]-1; i++) aj[i] -= 1; 1425 } 1426 for (i=0; i<ai[m]+shift; i++) col[aj[i]] += 1; 1427 ierr = MatCreateSeqAIJ(A->comm,A->n,m,0,col,&C);CHKERRQ(ierr); 1428 ierr = PetscFree(col);CHKERRQ(ierr); 1429 for (i=0; i<m; i++) { 1430 len = ai[i+1]-ai[i]; 1431 ierr = MatSetValues(C,len,aj,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 1432 array += len; 1433 aj += len; 1434 } 1435 if (shift) { 1436 for (i=0; i<ai[m]-1; i++) aj[i] += 1; 1437 } 1438 1439 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1440 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1441 1442 if (B) { 1443 *B = C; 1444 } else { 1445 ierr = MatHeaderCopy(A,C);CHKERRQ(ierr); 1446 } 1447 PetscFunctionReturn(0); 1448 } 1449 1450 #undef __FUNCT__ 1451 #define __FUNCT__ "MatDiagonalScale_SeqAIJ" 1452 int MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 1453 { 1454 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1455 PetscScalar *l,*r,x,*v; 1456 int ierr,i,j,m = A->m,n = A->n,M,nz = a->nz,*jj,shift = a->indexshift; 1457 1458 PetscFunctionBegin; 1459 if (ll) { 1460 /* The local size is used so that VecMPI can be passed to this routine 1461 by MatDiagonalScale_MPIAIJ */ 1462 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 1463 if (m != A->m) SETERRQ(PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 1464 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1465 v = a->a; 1466 for (i=0; i<m; i++) { 1467 x = l[i]; 1468 M = a->i[i+1] - a->i[i]; 1469 for (j=0; j<M; j++) { (*v++) *= x;} 1470 } 1471 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1472 PetscLogFlops(nz); 1473 } 1474 if (rr) { 1475 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 1476 if (n != A->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 1477 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1478 v = a->a; jj = a->j; 1479 for (i=0; i<nz; i++) { 1480 (*v++) *= r[*jj++ + shift]; 1481 } 1482 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1483 PetscLogFlops(nz); 1484 } 1485 PetscFunctionReturn(0); 1486 } 1487 1488 #undef __FUNCT__ 1489 #define __FUNCT__ "MatGetSubMatrix_SeqAIJ" 1490 int MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,int csize,MatReuse scall,Mat *B) 1491 { 1492 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 1493 int *smap,i,k,kstart,kend,ierr,oldcols = A->n,*lens; 1494 int row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 1495 int *irow,*icol,nrows,ncols,shift = a->indexshift,*ssmap; 1496 int *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 1497 PetscScalar *a_new,*mat_a; 1498 Mat C; 1499 PetscTruth stride; 1500 1501 PetscFunctionBegin; 1502 ierr = ISSorted(isrow,(PetscTruth*)&i);CHKERRQ(ierr); 1503 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"ISrow is not sorted"); 1504 ierr = ISSorted(iscol,(PetscTruth*)&i);CHKERRQ(ierr); 1505 if (!i) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"IScol is not sorted"); 1506 1507 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1508 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1509 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1510 1511 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 1512 ierr = ISStride(iscol,&stride);CHKERRQ(ierr); 1513 if (stride && step == 1) { 1514 /* special case of contiguous rows */ 1515 ierr = PetscMalloc((2*nrows+1)*sizeof(int),&lens);CHKERRQ(ierr); 1516 starts = lens + nrows; 1517 /* loop over new rows determining lens and starting points */ 1518 for (i=0; i<nrows; i++) { 1519 kstart = ai[irow[i]]+shift; 1520 kend = kstart + ailen[irow[i]]; 1521 for (k=kstart; k<kend; k++) { 1522 if (aj[k]+shift >= first) { 1523 starts[i] = k; 1524 break; 1525 } 1526 } 1527 sum = 0; 1528 while (k < kend) { 1529 if (aj[k++]+shift >= first+ncols) break; 1530 sum++; 1531 } 1532 lens[i] = sum; 1533 } 1534 /* create submatrix */ 1535 if (scall == MAT_REUSE_MATRIX) { 1536 int n_cols,n_rows; 1537 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1538 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 1539 ierr = MatZeroEntries(*B);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 1546 /* loop over rows inserting into submatrix */ 1547 a_new = c->a; 1548 j_new = c->j; 1549 i_new = c->i; 1550 i_new[0] = -shift; 1551 for (i=0; i<nrows; i++) { 1552 ii = starts[i]; 1553 lensi = lens[i]; 1554 for (k=0; k<lensi; k++) { 1555 *j_new++ = aj[ii+k] - first; 1556 } 1557 ierr = PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(PetscScalar));CHKERRQ(ierr); 1558 a_new += lensi; 1559 i_new[i+1] = i_new[i] + lensi; 1560 c->ilen[i] = lensi; 1561 } 1562 ierr = PetscFree(lens);CHKERRQ(ierr); 1563 } else { 1564 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1565 ierr = PetscMalloc((1+oldcols)*sizeof(int),&smap);CHKERRQ(ierr); 1566 ssmap = smap + shift; 1567 ierr = PetscMalloc((1+nrows)*sizeof(int),&lens);CHKERRQ(ierr); 1568 ierr = PetscMemzero(smap,oldcols*sizeof(int));CHKERRQ(ierr); 1569 for (i=0; i<ncols; i++) smap[icol[i]] = i+1; 1570 /* determine lens of each row */ 1571 for (i=0; i<nrows; i++) { 1572 kstart = ai[irow[i]]+shift; 1573 kend = kstart + a->ilen[irow[i]]; 1574 lens[i] = 0; 1575 for (k=kstart; k<kend; k++) { 1576 if (ssmap[aj[k]]) { 1577 lens[i]++; 1578 } 1579 } 1580 } 1581 /* Create and fill new matrix */ 1582 if (scall == MAT_REUSE_MATRIX) { 1583 PetscTruth equal; 1584 1585 c = (Mat_SeqAIJ *)((*B)->data); 1586 if ((*B)->m != nrows || (*B)->n != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 1587 ierr = PetscMemcmp(c->ilen,lens,(*B)->m*sizeof(int),&equal);CHKERRQ(ierr); 1588 if (!equal) { 1589 SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 1590 } 1591 ierr = PetscMemzero(c->ilen,(*B)->m*sizeof(int));CHKERRQ(ierr); 1592 C = *B; 1593 } else { 1594 ierr = MatCreateSeqAIJ(A->comm,nrows,ncols,0,lens,&C);CHKERRQ(ierr); 1595 } 1596 c = (Mat_SeqAIJ *)(C->data); 1597 for (i=0; i<nrows; i++) { 1598 row = irow[i]; 1599 kstart = ai[row]+shift; 1600 kend = kstart + a->ilen[row]; 1601 mat_i = c->i[i]+shift; 1602 mat_j = c->j + mat_i; 1603 mat_a = c->a + mat_i; 1604 mat_ilen = c->ilen + i; 1605 for (k=kstart; k<kend; k++) { 1606 if ((tcol=ssmap[a->j[k]])) { 1607 *mat_j++ = tcol - (!shift); 1608 *mat_a++ = a->a[k]; 1609 (*mat_ilen)++; 1610 1611 } 1612 } 1613 } 1614 /* Free work space */ 1615 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1616 ierr = PetscFree(smap);CHKERRQ(ierr); 1617 ierr = PetscFree(lens);CHKERRQ(ierr); 1618 } 1619 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1620 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1621 1622 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1623 *B = C; 1624 PetscFunctionReturn(0); 1625 } 1626 1627 /* 1628 */ 1629 #undef __FUNCT__ 1630 #define __FUNCT__ "MatILUFactor_SeqAIJ" 1631 int MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,MatILUInfo *info) 1632 { 1633 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1634 int ierr; 1635 Mat outA; 1636 PetscTruth row_identity,col_identity; 1637 1638 PetscFunctionBegin; 1639 if (info && info->levels != 0) SETERRQ(PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 1640 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 1641 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 1642 if (!row_identity || !col_identity) { 1643 SETERRQ(1,"Row and column permutations must be identity for in-place ILU"); 1644 } 1645 1646 outA = inA; 1647 inA->factor = FACTOR_LU; 1648 a->row = row; 1649 a->col = col; 1650 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1651 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 1652 1653 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 1654 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} /* need to remove old one */ 1655 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 1656 PetscLogObjectParent(inA,a->icol); 1657 1658 if (!a->solve_work) { /* this matrix may have been factored before */ 1659 ierr = PetscMalloc((inA->m+1)*sizeof(PetscScalar),&a->solve_work);CHKERRQ(ierr); 1660 } 1661 1662 if (!a->diag) { 1663 ierr = MatMarkDiagonal_SeqAIJ(inA);CHKERRQ(ierr); 1664 } 1665 ierr = MatLUFactorNumeric_SeqAIJ(inA,&outA);CHKERRQ(ierr); 1666 PetscFunctionReturn(0); 1667 } 1668 1669 #include "petscblaslapack.h" 1670 #undef __FUNCT__ 1671 #define __FUNCT__ "MatScale_SeqAIJ" 1672 int MatScale_SeqAIJ(PetscScalar *alpha,Mat inA) 1673 { 1674 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1675 int one = 1; 1676 1677 PetscFunctionBegin; 1678 BLscal_(&a->nz,alpha,a->a,&one); 1679 PetscLogFlops(a->nz); 1680 PetscFunctionReturn(0); 1681 } 1682 1683 #undef __FUNCT__ 1684 #define __FUNCT__ "MatGetSubMatrices_SeqAIJ" 1685 int MatGetSubMatrices_SeqAIJ(Mat A,int n,IS *irow,IS *icol,MatReuse scall,Mat **B) 1686 { 1687 int ierr,i; 1688 1689 PetscFunctionBegin; 1690 if (scall == MAT_INITIAL_MATRIX) { 1691 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1692 } 1693 1694 for (i=0; i<n; i++) { 1695 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1696 } 1697 PetscFunctionReturn(0); 1698 } 1699 1700 #undef __FUNCT__ 1701 #define __FUNCT__ "MatGetBlockSize_SeqAIJ" 1702 int MatGetBlockSize_SeqAIJ(Mat A,int *bs) 1703 { 1704 PetscFunctionBegin; 1705 *bs = 1; 1706 PetscFunctionReturn(0); 1707 } 1708 1709 #undef __FUNCT__ 1710 #define __FUNCT__ "MatIncreaseOverlap_SeqAIJ" 1711 int MatIncreaseOverlap_SeqAIJ(Mat A,int is_max,IS *is,int ov) 1712 { 1713 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1714 int shift,row,i,j,k,l,m,n,*idx,ierr,*nidx,isz,val; 1715 int start,end,*ai,*aj; 1716 PetscBT table; 1717 1718 PetscFunctionBegin; 1719 shift = a->indexshift; 1720 m = A->m; 1721 ai = a->i; 1722 aj = a->j+shift; 1723 1724 if (ov < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"illegal overlap value used"); 1725 1726 ierr = PetscMalloc((m+1)*sizeof(int),&nidx);CHKERRQ(ierr); 1727 ierr = PetscBTCreate(m,table);CHKERRQ(ierr); 1728 1729 for (i=0; i<is_max; i++) { 1730 /* Initialize the two local arrays */ 1731 isz = 0; 1732 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 1733 1734 /* Extract the indices, assume there can be duplicate entries */ 1735 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 1736 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 1737 1738 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 1739 for (j=0; j<n ; ++j){ 1740 if(!PetscBTLookupSet(table,idx[j])) { nidx[isz++] = idx[j];} 1741 } 1742 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 1743 ierr = ISDestroy(is[i]);CHKERRQ(ierr); 1744 1745 k = 0; 1746 for (j=0; j<ov; j++){ /* for each overlap */ 1747 n = isz; 1748 for (; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 1749 row = nidx[k]; 1750 start = ai[row]; 1751 end = ai[row+1]; 1752 for (l = start; l<end ; l++){ 1753 val = aj[l] + shift; 1754 if (!PetscBTLookupSet(table,val)) {nidx[isz++] = val;} 1755 } 1756 } 1757 } 1758 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,(is+i));CHKERRQ(ierr); 1759 } 1760 ierr = PetscBTDestroy(table);CHKERRQ(ierr); 1761 ierr = PetscFree(nidx);CHKERRQ(ierr); 1762 PetscFunctionReturn(0); 1763 } 1764 1765 /* -------------------------------------------------------------- */ 1766 #undef __FUNCT__ 1767 #define __FUNCT__ "MatPermute_SeqAIJ" 1768 int MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 1769 { 1770 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1771 PetscScalar *vwork; 1772 int i,ierr,nz,m = A->m,n = A->n,*cwork; 1773 int *row,*col,*cnew,j,*lens; 1774 IS icolp,irowp; 1775 1776 PetscFunctionBegin; 1777 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 1778 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 1779 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 1780 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 1781 1782 /* determine lengths of permuted rows */ 1783 ierr = PetscMalloc((m+1)*sizeof(int),&lens);CHKERRQ(ierr); 1784 for (i=0; i<m; i++) { 1785 lens[row[i]] = a->i[i+1] - a->i[i]; 1786 } 1787 ierr = MatCreateSeqAIJ(A->comm,m,n,0,lens,B);CHKERRQ(ierr); 1788 ierr = PetscFree(lens);CHKERRQ(ierr); 1789 1790 ierr = PetscMalloc(n*sizeof(int),&cnew);CHKERRQ(ierr); 1791 for (i=0; i<m; i++) { 1792 ierr = MatGetRow(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 1793 for (j=0; j<nz; j++) { cnew[j] = col[cwork[j]];} 1794 ierr = MatSetValues(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 1795 ierr = MatRestoreRow(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 1796 } 1797 ierr = PetscFree(cnew);CHKERRQ(ierr); 1798 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1799 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1800 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 1801 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 1802 ierr = ISDestroy(irowp);CHKERRQ(ierr); 1803 ierr = ISDestroy(icolp);CHKERRQ(ierr); 1804 PetscFunctionReturn(0); 1805 } 1806 1807 #undef __FUNCT__ 1808 #define __FUNCT__ "MatPrintHelp_SeqAIJ" 1809 int MatPrintHelp_SeqAIJ(Mat A) 1810 { 1811 static PetscTruth called = PETSC_FALSE; 1812 MPI_Comm comm = A->comm; 1813 int ierr; 1814 1815 PetscFunctionBegin; 1816 if (called) {PetscFunctionReturn(0);} else called = PETSC_TRUE; 1817 ierr = (*PetscHelpPrintf)(comm," Options for MATSEQAIJ and MATMPIAIJ matrix formats (the defaults):\n");CHKERRQ(ierr); 1818 ierr = (*PetscHelpPrintf)(comm," -mat_lu_pivotthreshold <threshold>: Set pivoting threshold\n");CHKERRQ(ierr); 1819 ierr = (*PetscHelpPrintf)(comm," -mat_aij_oneindex: internal indices begin at 1 instead of the default 0.\n");CHKERRQ(ierr); 1820 ierr = (*PetscHelpPrintf)(comm," -mat_aij_no_inode: Do not use inodes\n");CHKERRQ(ierr); 1821 ierr = (*PetscHelpPrintf)(comm," -mat_aij_inode_limit <limit>: Set inode limit (max limit=5)\n");CHKERRQ(ierr); 1822 #if defined(PETSC_HAVE_ESSL) 1823 ierr = (*PetscHelpPrintf)(comm," -mat_aij_essl: Use IBM sparse LU factorization and solve.\n");CHKERRQ(ierr); 1824 #endif 1825 #if defined(PETSC_HAVE_LUSOL) 1826 ierr = (*PetscHelpPrintf)(comm," -mat_aij_lusol: Use the Stanford LUSOL sparse factorization and solve.\n");CHKERRQ(ierr); 1827 #endif 1828 #if defined(PETSC_HAVE_MATLAB_ENGINE) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) 1829 ierr = (*PetscHelpPrintf)(comm," -mat_aij_matlab: Use Matlab engine sparse LU factorization and solve.\n");CHKERRQ(ierr); 1830 #endif 1831 PetscFunctionReturn(0); 1832 } 1833 EXTERN int MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg); 1834 EXTERN int MatFDColoringCreate_SeqAIJ(Mat,ISColoring,MatFDColoring); 1835 EXTERN int MatILUDTFactor_SeqAIJ(Mat,MatILUInfo*,IS,IS,Mat*); 1836 #undef __FUNCT__ 1837 #define __FUNCT__ "MatCopy_SeqAIJ" 1838 int MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 1839 { 1840 int ierr; 1841 PetscTruth flg; 1842 1843 PetscFunctionBegin; 1844 ierr = PetscTypeCompare((PetscObject)B,MATSEQAIJ,&flg);CHKERRQ(ierr); 1845 if (str == SAME_NONZERO_PATTERN && flg) { 1846 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1847 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 1848 1849 if (a->i[A->m]+a->indexshift != b->i[B->m]+a->indexshift) { 1850 SETERRQ(1,"Number of nonzeros in two matrices are different"); 1851 } 1852 ierr = PetscMemcpy(b->a,a->a,(a->i[A->m]+a->indexshift)*sizeof(PetscScalar));CHKERRQ(ierr); 1853 } else { 1854 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1855 } 1856 PetscFunctionReturn(0); 1857 } 1858 1859 #undef __FUNCT__ 1860 #define __FUNCT__ "MatSetUpPreallocation_SeqAIJ" 1861 int MatSetUpPreallocation_SeqAIJ(Mat A) 1862 { 1863 int ierr; 1864 1865 PetscFunctionBegin; 1866 ierr = MatSeqAIJSetPreallocation(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 1867 PetscFunctionReturn(0); 1868 } 1869 1870 #undef __FUNCT__ 1871 #define __FUNCT__ "MatGetArray_SeqAIJ" 1872 int MatGetArray_SeqAIJ(Mat A,PetscScalar **array) 1873 { 1874 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1875 PetscFunctionBegin; 1876 *array = a->a; 1877 PetscFunctionReturn(0); 1878 } 1879 1880 #undef __FUNCT__ 1881 #define __FUNCT__ "MatRestoreArray_SeqAIJ" 1882 int MatRestoreArray_SeqAIJ(Mat A,PetscScalar **array) 1883 { 1884 PetscFunctionBegin; 1885 PetscFunctionReturn(0); 1886 } 1887 1888 #undef __FUNCT__ 1889 #define __FUNCT__ "MatFDColoringApply_SeqAIJ" 1890 int MatFDColoringApply_SeqAIJ(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx) 1891 { 1892 int (*f)(void *,Vec,Vec,void*) = (int (*)(void *,Vec,Vec,void *))coloring->f; 1893 int k,ierr,N,start,end,l,row,col,srow,**vscaleforrow,m1,m2; 1894 PetscScalar dx,mone = -1.0,*y,*xx,*w3_array; 1895 PetscScalar *vscale_array; 1896 PetscReal epsilon = coloring->error_rel,umin = coloring->umin; 1897 Vec w1,w2,w3; 1898 void *fctx = coloring->fctx; 1899 PetscTruth flg; 1900 1901 PetscFunctionBegin; 1902 if (!coloring->w1) { 1903 ierr = VecDuplicate(x1,&coloring->w1);CHKERRQ(ierr); 1904 PetscLogObjectParent(coloring,coloring->w1); 1905 ierr = VecDuplicate(x1,&coloring->w2);CHKERRQ(ierr); 1906 PetscLogObjectParent(coloring,coloring->w2); 1907 ierr = VecDuplicate(x1,&coloring->w3);CHKERRQ(ierr); 1908 PetscLogObjectParent(coloring,coloring->w3); 1909 } 1910 w1 = coloring->w1; w2 = coloring->w2; w3 = coloring->w3; 1911 1912 ierr = MatSetUnfactored(J);CHKERRQ(ierr); 1913 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_fd_coloring_dont_rezero",&flg);CHKERRQ(ierr); 1914 if (flg) { 1915 PetscLogInfo(coloring,"MatFDColoringApply_SeqAIJ: Not calling MatZeroEntries()\n"); 1916 } else { 1917 ierr = MatZeroEntries(J);CHKERRQ(ierr); 1918 } 1919 1920 ierr = VecGetOwnershipRange(x1,&start,&end);CHKERRQ(ierr); 1921 ierr = VecGetSize(x1,&N);CHKERRQ(ierr); 1922 1923 /* 1924 This is a horrible, horrible, hack. See DMMGComputeJacobian_Multigrid() it inproperly sets 1925 coloring->F for the coarser grids from the finest 1926 */ 1927 if (coloring->F) { 1928 ierr = VecGetLocalSize(coloring->F,&m1);CHKERRQ(ierr); 1929 ierr = VecGetLocalSize(w1,&m2);CHKERRQ(ierr); 1930 if (m1 != m2) { 1931 coloring->F = 0; 1932 } 1933 } 1934 1935 if (coloring->F) { 1936 w1 = coloring->F; 1937 coloring->F = 0; 1938 } else { 1939 ierr = (*f)(sctx,x1,w1,fctx);CHKERRQ(ierr); 1940 } 1941 1942 /* 1943 Compute all the scale factors and share with other processors 1944 */ 1945 ierr = VecGetArray(x1,&xx);CHKERRQ(ierr);xx = xx - start; 1946 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr);vscale_array = vscale_array - start; 1947 for (k=0; k<coloring->ncolors; k++) { 1948 /* 1949 Loop over each column associated with color adding the 1950 perturbation to the vector w3. 1951 */ 1952 for (l=0; l<coloring->ncolumns[k]; l++) { 1953 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 1954 dx = xx[col]; 1955 if (dx == 0.0) dx = 1.0; 1956 #if !defined(PETSC_USE_COMPLEX) 1957 if (dx < umin && dx >= 0.0) dx = umin; 1958 else if (dx < 0.0 && dx > -umin) dx = -umin; 1959 #else 1960 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 1961 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 1962 #endif 1963 dx *= epsilon; 1964 vscale_array[col] = 1.0/dx; 1965 } 1966 } 1967 vscale_array = vscale_array + start;ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 1968 ierr = VecGhostUpdateBegin(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1969 ierr = VecGhostUpdateEnd(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1970 1971 /* ierr = VecView(coloring->vscale,PETSC_VIEWER_STDOUT_WORLD); 1972 ierr = VecView(x1,PETSC_VIEWER_STDOUT_WORLD);*/ 1973 1974 if (coloring->vscaleforrow) vscaleforrow = coloring->vscaleforrow; 1975 else vscaleforrow = coloring->columnsforrow; 1976 1977 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 1978 /* 1979 Loop over each color 1980 */ 1981 for (k=0; k<coloring->ncolors; k++) { 1982 ierr = VecCopy(x1,w3);CHKERRQ(ierr); 1983 ierr = VecGetArray(w3,&w3_array);CHKERRQ(ierr);w3_array = w3_array - start; 1984 /* 1985 Loop over each column associated with color adding the 1986 perturbation to the vector w3. 1987 */ 1988 for (l=0; l<coloring->ncolumns[k]; l++) { 1989 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 1990 dx = xx[col]; 1991 if (dx == 0.0) dx = 1.0; 1992 #if !defined(PETSC_USE_COMPLEX) 1993 if (dx < umin && dx >= 0.0) dx = umin; 1994 else if (dx < 0.0 && dx > -umin) dx = -umin; 1995 #else 1996 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 1997 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 1998 #endif 1999 dx *= epsilon; 2000 if (!PetscAbsScalar(dx)) SETERRQ(1,"Computed 0 differencing parameter"); 2001 w3_array[col] += dx; 2002 } 2003 w3_array = w3_array + start; ierr = VecRestoreArray(w3,&w3_array);CHKERRQ(ierr); 2004 2005 /* 2006 Evaluate function at x1 + dx (here dx is a vector of perturbations) 2007 */ 2008 2009 ierr = (*f)(sctx,w3,w2,fctx);CHKERRQ(ierr); 2010 ierr = VecAXPY(&mone,w1,w2);CHKERRQ(ierr); 2011 2012 /* 2013 Loop over rows of vector, putting results into Jacobian matrix 2014 */ 2015 ierr = VecGetArray(w2,&y);CHKERRQ(ierr); 2016 for (l=0; l<coloring->nrows[k]; l++) { 2017 row = coloring->rows[k][l]; 2018 col = coloring->columnsforrow[k][l]; 2019 y[row] *= vscale_array[vscaleforrow[k][l]]; 2020 srow = row + start; 2021 ierr = MatSetValues_SeqAIJ(J,1,&srow,1,&col,y+row,INSERT_VALUES);CHKERRQ(ierr); 2022 } 2023 ierr = VecRestoreArray(w2,&y);CHKERRQ(ierr); 2024 } 2025 ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2026 xx = xx + start; ierr = VecRestoreArray(x1,&xx);CHKERRQ(ierr); 2027 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2028 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2029 PetscFunctionReturn(0); 2030 } 2031 2032 /* -------------------------------------------------------------------*/ 2033 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 2034 MatGetRow_SeqAIJ, 2035 MatRestoreRow_SeqAIJ, 2036 MatMult_SeqAIJ, 2037 MatMultAdd_SeqAIJ, 2038 MatMultTranspose_SeqAIJ, 2039 MatMultTransposeAdd_SeqAIJ, 2040 MatSolve_SeqAIJ, 2041 MatSolveAdd_SeqAIJ, 2042 MatSolveTranspose_SeqAIJ, 2043 MatSolveTransposeAdd_SeqAIJ, 2044 MatLUFactor_SeqAIJ, 2045 0, 2046 MatRelax_SeqAIJ, 2047 MatTranspose_SeqAIJ, 2048 MatGetInfo_SeqAIJ, 2049 MatEqual_SeqAIJ, 2050 MatGetDiagonal_SeqAIJ, 2051 MatDiagonalScale_SeqAIJ, 2052 MatNorm_SeqAIJ, 2053 0, 2054 MatAssemblyEnd_SeqAIJ, 2055 MatCompress_SeqAIJ, 2056 MatSetOption_SeqAIJ, 2057 MatZeroEntries_SeqAIJ, 2058 MatZeroRows_SeqAIJ, 2059 MatLUFactorSymbolic_SeqAIJ, 2060 MatLUFactorNumeric_SeqAIJ, 2061 0, 2062 0, 2063 MatSetUpPreallocation_SeqAIJ, 2064 0, 2065 MatGetOwnershipRange_SeqAIJ, 2066 MatILUFactorSymbolic_SeqAIJ, 2067 0, 2068 MatGetArray_SeqAIJ, 2069 MatRestoreArray_SeqAIJ, 2070 MatDuplicate_SeqAIJ, 2071 0, 2072 0, 2073 MatILUFactor_SeqAIJ, 2074 0, 2075 0, 2076 MatGetSubMatrices_SeqAIJ, 2077 MatIncreaseOverlap_SeqAIJ, 2078 MatGetValues_SeqAIJ, 2079 MatCopy_SeqAIJ, 2080 MatPrintHelp_SeqAIJ, 2081 MatScale_SeqAIJ, 2082 0, 2083 0, 2084 MatILUDTFactor_SeqAIJ, 2085 MatGetBlockSize_SeqAIJ, 2086 MatGetRowIJ_SeqAIJ, 2087 MatRestoreRowIJ_SeqAIJ, 2088 MatGetColumnIJ_SeqAIJ, 2089 MatRestoreColumnIJ_SeqAIJ, 2090 MatFDColoringCreate_SeqAIJ, 2091 0, 2092 0, 2093 MatPermute_SeqAIJ, 2094 0, 2095 0, 2096 MatDestroy_SeqAIJ, 2097 MatView_SeqAIJ, 2098 MatGetPetscMaps_Petsc, 2099 0, 2100 0, 2101 0, 2102 0, 2103 0, 2104 0, 2105 0, 2106 0, 2107 MatSetColoring_SeqAIJ, 2108 MatSetValuesAdic_SeqAIJ, 2109 MatSetValuesAdifor_SeqAIJ, 2110 MatFDColoringApply_SeqAIJ}; 2111 2112 EXTERN int MatUseSuperLU_SeqAIJ(Mat); 2113 EXTERN int MatUseEssl_SeqAIJ(Mat); 2114 EXTERN int MatUseLUSOL_SeqAIJ(Mat); 2115 EXTERN int MatUseMatlab_SeqAIJ(Mat); 2116 EXTERN int MatUseDXML_SeqAIJ(Mat); 2117 2118 EXTERN_C_BEGIN 2119 #undef __FUNCT__ 2120 #define __FUNCT__ "MatSeqAIJSetColumnIndices_SeqAIJ" 2121 2122 int MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,int *indices) 2123 { 2124 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2125 int i,nz,n; 2126 2127 PetscFunctionBegin; 2128 2129 nz = aij->maxnz; 2130 n = mat->n; 2131 for (i=0; i<nz; i++) { 2132 aij->j[i] = indices[i]; 2133 } 2134 aij->nz = nz; 2135 for (i=0; i<n; i++) { 2136 aij->ilen[i] = aij->imax[i]; 2137 } 2138 2139 PetscFunctionReturn(0); 2140 } 2141 EXTERN_C_END 2142 2143 #undef __FUNCT__ 2144 #define __FUNCT__ "MatSeqAIJSetColumnIndices" 2145 /*@ 2146 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 2147 in the matrix. 2148 2149 Input Parameters: 2150 + mat - the SeqAIJ matrix 2151 - indices - the column indices 2152 2153 Level: advanced 2154 2155 Notes: 2156 This can be called if you have precomputed the nonzero structure of the 2157 matrix and want to provide it to the matrix object to improve the performance 2158 of the MatSetValues() operation. 2159 2160 You MUST have set the correct numbers of nonzeros per row in the call to 2161 MatCreateSeqAIJ(). 2162 2163 MUST be called before any calls to MatSetValues(); 2164 2165 The indices should start with zero, not one. 2166 2167 @*/ 2168 int MatSeqAIJSetColumnIndices(Mat mat,int *indices) 2169 { 2170 int ierr,(*f)(Mat,int *); 2171 2172 PetscFunctionBegin; 2173 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2174 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqAIJSetColumnIndices_C",(void (**)())&f);CHKERRQ(ierr); 2175 if (f) { 2176 ierr = (*f)(mat,indices);CHKERRQ(ierr); 2177 } else { 2178 SETERRQ(1,"Wrong type of matrix to set column indices"); 2179 } 2180 PetscFunctionReturn(0); 2181 } 2182 2183 /* ----------------------------------------------------------------------------------------*/ 2184 2185 EXTERN_C_BEGIN 2186 #undef __FUNCT__ 2187 #define __FUNCT__ "MatStoreValues_SeqAIJ" 2188 int MatStoreValues_SeqAIJ(Mat mat) 2189 { 2190 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2191 int nz = aij->i[mat->m]+aij->indexshift,ierr; 2192 2193 PetscFunctionBegin; 2194 if (aij->nonew != 1) { 2195 SETERRQ(1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 2196 } 2197 2198 /* allocate space for values if not already there */ 2199 if (!aij->saved_values) { 2200 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 2201 } 2202 2203 /* copy values over */ 2204 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2205 PetscFunctionReturn(0); 2206 } 2207 EXTERN_C_END 2208 2209 #undef __FUNCT__ 2210 #define __FUNCT__ "MatStoreValues" 2211 /*@ 2212 MatStoreValues - Stashes a copy of the matrix values; this allows, for 2213 example, reuse of the linear part of a Jacobian, while recomputing the 2214 nonlinear portion. 2215 2216 Collect on Mat 2217 2218 Input Parameters: 2219 . mat - the matrix (currently on AIJ matrices support this option) 2220 2221 Level: advanced 2222 2223 Common Usage, with SNESSolve(): 2224 $ Create Jacobian matrix 2225 $ Set linear terms into matrix 2226 $ Apply boundary conditions to matrix, at this time matrix must have 2227 $ final nonzero structure (i.e. setting the nonlinear terms and applying 2228 $ boundary conditions again will not change the nonzero structure 2229 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 2230 $ ierr = MatStoreValues(mat); 2231 $ Call SNESSetJacobian() with matrix 2232 $ In your Jacobian routine 2233 $ ierr = MatRetrieveValues(mat); 2234 $ Set nonlinear terms in matrix 2235 2236 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 2237 $ // build linear portion of Jacobian 2238 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 2239 $ ierr = MatStoreValues(mat); 2240 $ loop over nonlinear iterations 2241 $ ierr = MatRetrieveValues(mat); 2242 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 2243 $ // call MatAssemblyBegin/End() on matrix 2244 $ Solve linear system with Jacobian 2245 $ endloop 2246 2247 Notes: 2248 Matrix must already be assemblied before calling this routine 2249 Must set the matrix option MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); before 2250 calling this routine. 2251 2252 .seealso: MatRetrieveValues() 2253 2254 @*/ 2255 int MatStoreValues(Mat mat) 2256 { 2257 int ierr,(*f)(Mat); 2258 2259 PetscFunctionBegin; 2260 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2261 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2262 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2263 2264 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatStoreValues_C",(void (**)())&f);CHKERRQ(ierr); 2265 if (f) { 2266 ierr = (*f)(mat);CHKERRQ(ierr); 2267 } else { 2268 SETERRQ(1,"Wrong type of matrix to store values"); 2269 } 2270 PetscFunctionReturn(0); 2271 } 2272 2273 EXTERN_C_BEGIN 2274 #undef __FUNCT__ 2275 #define __FUNCT__ "MatRetrieveValues_SeqAIJ" 2276 int MatRetrieveValues_SeqAIJ(Mat mat) 2277 { 2278 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2279 int nz = aij->i[mat->m]+aij->indexshift,ierr; 2280 2281 PetscFunctionBegin; 2282 if (aij->nonew != 1) { 2283 SETERRQ(1,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 2284 } 2285 if (!aij->saved_values) { 2286 SETERRQ(1,"Must call MatStoreValues(A);first"); 2287 } 2288 2289 /* copy values over */ 2290 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2291 PetscFunctionReturn(0); 2292 } 2293 EXTERN_C_END 2294 2295 #undef __FUNCT__ 2296 #define __FUNCT__ "MatRetrieveValues" 2297 /*@ 2298 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 2299 example, reuse of the linear part of a Jacobian, while recomputing the 2300 nonlinear portion. 2301 2302 Collect on Mat 2303 2304 Input Parameters: 2305 . mat - the matrix (currently on AIJ matrices support this option) 2306 2307 Level: advanced 2308 2309 .seealso: MatStoreValues() 2310 2311 @*/ 2312 int MatRetrieveValues(Mat mat) 2313 { 2314 int ierr,(*f)(Mat); 2315 2316 PetscFunctionBegin; 2317 PetscValidHeaderSpecific(mat,MAT_COOKIE); 2318 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2319 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2320 2321 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatRetrieveValues_C",(void (**)())&f);CHKERRQ(ierr); 2322 if (f) { 2323 ierr = (*f)(mat);CHKERRQ(ierr); 2324 } else { 2325 SETERRQ(1,"Wrong type of matrix to retrieve values"); 2326 } 2327 PetscFunctionReturn(0); 2328 } 2329 2330 /* 2331 This allows SeqAIJ matrices to be passed to the matlab engine 2332 */ 2333 #if defined(PETSC_HAVE_MATLAB_ENGINE) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) 2334 #include "engine.h" /* Matlab include file */ 2335 #include "mex.h" /* Matlab include file */ 2336 EXTERN_C_BEGIN 2337 #undef __FUNCT__ 2338 #define __FUNCT__ "MatMatlabEnginePut_SeqAIJ" 2339 int MatMatlabEnginePut_SeqAIJ(PetscObject obj,void *engine) 2340 { 2341 int ierr,i,*ai,*aj; 2342 Mat B = (Mat)obj; 2343 PetscScalar *array; 2344 mxArray *mat; 2345 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)B->data; 2346 2347 PetscFunctionBegin; 2348 mat = mxCreateSparse(B->n,B->m,aij->nz,mxREAL); 2349 ierr = PetscMemcpy(mxGetPr(mat),aij->a,aij->nz*sizeof(PetscScalar));CHKERRQ(ierr); 2350 /* Matlab stores by column, not row so we pass in the transpose of the matrix */ 2351 ierr = PetscMemcpy(mxGetIr(mat),aij->j,aij->nz*sizeof(int));CHKERRQ(ierr); 2352 ierr = PetscMemcpy(mxGetJc(mat),aij->i,(B->m+1)*sizeof(int));CHKERRQ(ierr); 2353 2354 /* Matlab indices start at 0 for sparse (what a surprise) */ 2355 if (aij->indexshift) { 2356 for (i=0; i<B->m+1; i++) { 2357 ai[i]--; 2358 } 2359 for (i=0; i<aij->nz; i++) { 2360 aj[i]--; 2361 } 2362 } 2363 ierr = PetscObjectName(obj);CHKERRQ(ierr); 2364 mxSetName(mat,obj->name); 2365 engPutArray((Engine *)engine,mat); 2366 PetscFunctionReturn(0); 2367 } 2368 EXTERN_C_END 2369 2370 EXTERN_C_BEGIN 2371 #undef __FUNCT__ 2372 #define __FUNCT__ "MatMatlabEngineGet_SeqAIJ" 2373 int MatMatlabEngineGet_SeqAIJ(PetscObject obj,void *engine) 2374 { 2375 int ierr,ii; 2376 Mat mat = (Mat)obj; 2377 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 2378 mxArray *mmat; 2379 2380 PetscFunctionBegin; 2381 ierr = PetscFree(aij->a);CHKERRQ(ierr); 2382 aij->indexshift = 0; 2383 2384 mmat = engGetArray((Engine *)engine,obj->name); 2385 2386 aij->nz = (mxGetJc(mmat))[mat->m]; 2387 ierr = PetscMalloc(aij->nz*(sizeof(int)+sizeof(PetscScalar))+(mat->m+1)*sizeof(int),&aij->a);CHKERRQ(ierr); 2388 aij->j = (int*)(aij->a + aij->nz); 2389 aij->i = aij->j + aij->nz; 2390 aij->singlemalloc = PETSC_TRUE; 2391 aij->freedata = PETSC_TRUE; 2392 2393 ierr = PetscMemcpy(aij->a,mxGetPr(mmat),aij->nz*sizeof(PetscScalar));CHKERRQ(ierr); 2394 /* Matlab stores by column, not row so we pass in the transpose of the matrix */ 2395 ierr = PetscMemcpy(aij->j,mxGetIr(mmat),aij->nz*sizeof(int));CHKERRQ(ierr); 2396 ierr = PetscMemcpy(aij->i,mxGetJc(mmat),(mat->m+1)*sizeof(int));CHKERRQ(ierr); 2397 2398 for (ii=0; ii<mat->m; ii++) { 2399 aij->ilen[ii] = aij->imax[ii] = aij->i[ii+1] - aij->i[ii]; 2400 } 2401 2402 ierr = MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2403 ierr = MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2404 2405 PetscFunctionReturn(0); 2406 } 2407 EXTERN_C_END 2408 #endif 2409 2410 /* --------------------------------------------------------------------------------*/ 2411 2412 #undef __FUNCT__ 2413 #define __FUNCT__ "MatCreateSeqAIJ" 2414 /*@C 2415 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 2416 (the default parallel PETSc format). For good matrix assembly performance 2417 the user should preallocate the matrix storage by setting the parameter nz 2418 (or the array nnz). By setting these parameters accurately, performance 2419 during matrix assembly can be increased by more than a factor of 50. 2420 2421 Collective on MPI_Comm 2422 2423 Input Parameters: 2424 + comm - MPI communicator, set to PETSC_COMM_SELF 2425 . m - number of rows 2426 . n - number of columns 2427 . nz - number of nonzeros per row (same for all rows) 2428 - nnz - array containing the number of nonzeros in the various rows 2429 (possibly different for each row) or PETSC_NULL 2430 2431 Output Parameter: 2432 . A - the matrix 2433 2434 Notes: 2435 The AIJ format (also called the Yale sparse matrix format or 2436 compressed row storage), is fully compatible with standard Fortran 77 2437 storage. That is, the stored row and column indices can begin at 2438 either one (as in Fortran) or zero. See the users' manual for details. 2439 2440 Specify the preallocated storage with either nz or nnz (not both). 2441 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2442 allocation. For large problems you MUST preallocate memory or you 2443 will get TERRIBLE performance, see the users' manual chapter on matrices. 2444 2445 By default, this format uses inodes (identical nodes) when possible, to 2446 improve numerical efficiency of matrix-vector products and solves. We 2447 search for consecutive rows with the same nonzero structure, thereby 2448 reusing matrix information to achieve increased efficiency. 2449 2450 Options Database Keys: 2451 + -mat_aij_no_inode - Do not use inodes 2452 . -mat_aij_inode_limit <limit> - Sets inode limit (max limit=5) 2453 - -mat_aij_oneindex - Internally use indexing starting at 1 2454 rather than 0. Note that when calling MatSetValues(), 2455 the user still MUST index entries starting at 0! 2456 2457 Level: intermediate 2458 2459 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2460 2461 @*/ 2462 int MatCreateSeqAIJ(MPI_Comm comm,int m,int n,int nz,int *nnz,Mat *A) 2463 { 2464 int ierr; 2465 2466 PetscFunctionBegin; 2467 ierr = MatCreate(comm,m,n,m,n,A);CHKERRQ(ierr); 2468 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 2469 ierr = MatSeqAIJSetPreallocation(*A,nz,nnz);CHKERRQ(ierr); 2470 PetscFunctionReturn(0); 2471 } 2472 2473 #undef __FUNCT__ 2474 #define __FUNCT__ "MatSeqAIJSetPreallocation" 2475 /*@C 2476 MatSeqAIJSetPreallocation - For good matrix assembly performance 2477 the user should preallocate the matrix storage by setting the parameter nz 2478 (or the array nnz). By setting these parameters accurately, performance 2479 during matrix assembly can be increased by more than a factor of 50. 2480 2481 Collective on MPI_Comm 2482 2483 Input Parameters: 2484 + comm - MPI communicator, set to PETSC_COMM_SELF 2485 . m - number of rows 2486 . n - number of columns 2487 . nz - number of nonzeros per row (same for all rows) 2488 - nnz - array containing the number of nonzeros in the various rows 2489 (possibly different for each row) or PETSC_NULL 2490 2491 Output Parameter: 2492 . A - the matrix 2493 2494 Notes: 2495 The AIJ format (also called the Yale sparse matrix format or 2496 compressed row storage), is fully compatible with standard Fortran 77 2497 storage. That is, the stored row and column indices can begin at 2498 either one (as in Fortran) or zero. See the users' manual for details. 2499 2500 Specify the preallocated storage with either nz or nnz (not both). 2501 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2502 allocation. For large problems you MUST preallocate memory or you 2503 will get TERRIBLE performance, see the users' manual chapter on matrices. 2504 2505 By default, this format uses inodes (identical nodes) when possible, to 2506 improve numerical efficiency of matrix-vector products and solves. We 2507 search for consecutive rows with the same nonzero structure, thereby 2508 reusing matrix information to achieve increased efficiency. 2509 2510 Options Database Keys: 2511 + -mat_aij_no_inode - Do not use inodes 2512 . -mat_aij_inode_limit <limit> - Sets inode limit (max limit=5) 2513 - -mat_aij_oneindex - Internally use indexing starting at 1 2514 rather than 0. Note that when calling MatSetValues(), 2515 the user still MUST index entries starting at 0! 2516 2517 Level: intermediate 2518 2519 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2520 2521 @*/ 2522 int MatSeqAIJSetPreallocation(Mat B,int nz,int *nnz) 2523 { 2524 Mat_SeqAIJ *b; 2525 int i,len,ierr; 2526 PetscTruth flg2; 2527 2528 PetscFunctionBegin; 2529 ierr = PetscTypeCompare((PetscObject)B,MATSEQAIJ,&flg2);CHKERRQ(ierr); 2530 if (!flg2) PetscFunctionReturn(0); 2531 2532 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 2533 if (nz < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 2534 if (nnz) { 2535 for (i=0; i<B->m; i++) { 2536 if (nnz[i] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]); 2537 if (nnz[i] > B->n) SETERRQ3(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be greater than row length: local row %d value %d rowlength %d",i,nnz[i],B->n); 2538 } 2539 } 2540 2541 B->preallocated = PETSC_TRUE; 2542 b = (Mat_SeqAIJ*)B->data; 2543 2544 ierr = PetscMalloc((B->m+1)*sizeof(int),&b->imax);CHKERRQ(ierr); 2545 if (!nnz) { 2546 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 2547 else if (nz <= 0) nz = 1; 2548 for (i=0; i<B->m; i++) b->imax[i] = nz; 2549 nz = nz*B->m; 2550 } else { 2551 nz = 0; 2552 for (i=0; i<B->m; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 2553 } 2554 2555 /* allocate the matrix space */ 2556 len = nz*(sizeof(int) + sizeof(PetscScalar)) + (B->m+1)*sizeof(int); 2557 ierr = PetscMalloc(len,&b->a);CHKERRQ(ierr); 2558 b->j = (int*)(b->a + nz); 2559 ierr = PetscMemzero(b->j,nz*sizeof(int));CHKERRQ(ierr); 2560 b->i = b->j + nz; 2561 b->singlemalloc = PETSC_TRUE; 2562 b->freedata = PETSC_TRUE; 2563 2564 b->i[0] = -b->indexshift; 2565 for (i=1; i<B->m+1; i++) { 2566 b->i[i] = b->i[i-1] + b->imax[i-1]; 2567 } 2568 2569 /* b->ilen will count nonzeros in each row so far. */ 2570 ierr = PetscMalloc((B->m+1)*sizeof(int),&b->ilen);CHKERRQ(ierr); 2571 PetscLogObjectMemory(B,len+2*(B->m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2572 for (i=0; i<B->m; i++) { b->ilen[i] = 0;} 2573 2574 b->nz = 0; 2575 b->maxnz = nz; 2576 B->info.nz_unneeded = (double)b->maxnz; 2577 PetscFunctionReturn(0); 2578 } 2579 2580 EXTERN_C_BEGIN 2581 #undef __FUNCT__ 2582 #define __FUNCT__ "MatCreate_SeqAIJ" 2583 int MatCreate_SeqAIJ(Mat B) 2584 { 2585 Mat_SeqAIJ *b; 2586 int ierr,size; 2587 PetscTruth flg; 2588 2589 PetscFunctionBegin; 2590 ierr = MPI_Comm_size(B->comm,&size);CHKERRQ(ierr); 2591 if (size > 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 2592 2593 B->m = B->M = PetscMax(B->m,B->M); 2594 B->n = B->N = PetscMax(B->n,B->N); 2595 2596 ierr = PetscNew(Mat_SeqAIJ,&b);CHKERRQ(ierr); 2597 B->data = (void*)b; 2598 ierr = PetscMemzero(b,sizeof(Mat_SeqAIJ));CHKERRQ(ierr); 2599 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2600 B->factor = 0; 2601 B->lupivotthreshold = 1.0; 2602 B->mapping = 0; 2603 ierr = PetscOptionsGetReal(PETSC_NULL,"-mat_lu_pivotthreshold",&B->lupivotthreshold,PETSC_NULL);CHKERRQ(ierr); 2604 ierr = PetscOptionsHasName(PETSC_NULL,"-pc_ilu_preserve_row_sums",&b->ilu_preserve_row_sums);CHKERRQ(ierr); 2605 b->row = 0; 2606 b->col = 0; 2607 b->icol = 0; 2608 b->indexshift = 0; 2609 b->reallocs = 0; 2610 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_aij_oneindex",&flg);CHKERRQ(ierr); 2611 if (flg) b->indexshift = -1; 2612 2613 ierr = PetscMapCreateMPI(B->comm,B->m,B->m,&B->rmap);CHKERRQ(ierr); 2614 ierr = PetscMapCreateMPI(B->comm,B->n,B->n,&B->cmap);CHKERRQ(ierr); 2615 2616 b->sorted = PETSC_FALSE; 2617 b->ignorezeroentries = PETSC_FALSE; 2618 b->roworiented = PETSC_TRUE; 2619 b->nonew = 0; 2620 b->diag = 0; 2621 b->solve_work = 0; 2622 b->spptr = 0; 2623 b->inode.use = PETSC_TRUE; 2624 b->inode.node_count = 0; 2625 b->inode.size = 0; 2626 b->inode.limit = 5; 2627 b->inode.max_limit = 5; 2628 b->saved_values = 0; 2629 b->idiag = 0; 2630 b->ssor = 0; 2631 b->keepzeroedrows = PETSC_FALSE; 2632 2633 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 2634 2635 #if defined(PETSC_HAVE_SUPERLU) 2636 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_aij_superlu",&flg);CHKERRQ(ierr); 2637 if (flg) { ierr = MatUseSuperLU_SeqAIJ(B);CHKERRQ(ierr); } 2638 #endif 2639 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_aij_essl",&flg);CHKERRQ(ierr); 2640 if (flg) { ierr = MatUseEssl_SeqAIJ(B);CHKERRQ(ierr); } 2641 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_aij_lusol",&flg);CHKERRQ(ierr); 2642 if (flg) { ierr = MatUseLUSOL_SeqAIJ(B);CHKERRQ(ierr); } 2643 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_aij_matlab",&flg);CHKERRQ(ierr); 2644 if (flg) {ierr = MatUseMatlab_SeqAIJ(B);CHKERRQ(ierr);} 2645 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_aij_dxml",&flg);CHKERRQ(ierr); 2646 if (flg) { 2647 if (!b->indexshift) SETERRQ(PETSC_ERR_LIB,"need -mat_aij_oneindex with -mat_aij_dxml"); 2648 ierr = MatUseDXML_SeqAIJ(B);CHKERRQ(ierr); 2649 } 2650 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 2651 "MatSeqAIJSetColumnIndices_SeqAIJ", 2652 MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 2653 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 2654 "MatStoreValues_SeqAIJ", 2655 MatStoreValues_SeqAIJ);CHKERRQ(ierr); 2656 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 2657 "MatRetrieveValues_SeqAIJ", 2658 MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 2659 #if defined(PETSC_HAVE_MATLAB_ENGINE) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) 2660 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEnginePut_C","MatMatlabEnginePut_SeqAIJ",MatMatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 2661 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEngineGet_C","MatMatlabEngineGet_SeqAIJ",MatMatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 2662 #endif 2663 PetscFunctionReturn(0); 2664 } 2665 EXTERN_C_END 2666 2667 #undef __FUNCT__ 2668 #define __FUNCT__ "MatDuplicate_SeqAIJ" 2669 int MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 2670 { 2671 Mat C; 2672 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 2673 int i,len,m = A->m,shift = a->indexshift,ierr; 2674 2675 PetscFunctionBegin; 2676 *B = 0; 2677 ierr = MatCreate(A->comm,A->m,A->n,A->m,A->n,&C);CHKERRQ(ierr); 2678 ierr = MatSetType(C,MATSEQAIJ);CHKERRQ(ierr); 2679 c = (Mat_SeqAIJ*)C->data; 2680 2681 C->factor = A->factor; 2682 c->row = 0; 2683 c->col = 0; 2684 c->icol = 0; 2685 c->indexshift = shift; 2686 c->keepzeroedrows = a->keepzeroedrows; 2687 C->assembled = PETSC_TRUE; 2688 2689 C->M = A->m; 2690 C->N = A->n; 2691 2692 ierr = PetscMalloc((m+1)*sizeof(int),&c->imax);CHKERRQ(ierr); 2693 ierr = PetscMalloc((m+1)*sizeof(int),&c->ilen);CHKERRQ(ierr); 2694 for (i=0; i<m; i++) { 2695 c->imax[i] = a->imax[i]; 2696 c->ilen[i] = a->ilen[i]; 2697 } 2698 2699 /* allocate the matrix space */ 2700 c->singlemalloc = PETSC_TRUE; 2701 len = (m+1)*sizeof(int)+(a->i[m])*(sizeof(PetscScalar)+sizeof(int)); 2702 ierr = PetscMalloc(len,&c->a);CHKERRQ(ierr); 2703 c->j = (int*)(c->a + a->i[m] + shift); 2704 c->i = c->j + a->i[m] + shift; 2705 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(int));CHKERRQ(ierr); 2706 if (m > 0) { 2707 ierr = PetscMemcpy(c->j,a->j,(a->i[m]+shift)*sizeof(int));CHKERRQ(ierr); 2708 if (cpvalues == MAT_COPY_VALUES) { 2709 ierr = PetscMemcpy(c->a,a->a,(a->i[m]+shift)*sizeof(PetscScalar));CHKERRQ(ierr); 2710 } else { 2711 ierr = PetscMemzero(c->a,(a->i[m]+shift)*sizeof(PetscScalar));CHKERRQ(ierr); 2712 } 2713 } 2714 2715 PetscLogObjectMemory(C,len+2*(m+1)*sizeof(int)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2716 c->sorted = a->sorted; 2717 c->roworiented = a->roworiented; 2718 c->nonew = a->nonew; 2719 c->ilu_preserve_row_sums = a->ilu_preserve_row_sums; 2720 c->saved_values = 0; 2721 c->idiag = 0; 2722 c->ssor = 0; 2723 c->ignorezeroentries = a->ignorezeroentries; 2724 c->freedata = PETSC_TRUE; 2725 2726 if (a->diag) { 2727 ierr = PetscMalloc((m+1)*sizeof(int),&c->diag);CHKERRQ(ierr); 2728 PetscLogObjectMemory(C,(m+1)*sizeof(int)); 2729 for (i=0; i<m; i++) { 2730 c->diag[i] = a->diag[i]; 2731 } 2732 } else c->diag = 0; 2733 c->inode.use = a->inode.use; 2734 c->inode.limit = a->inode.limit; 2735 c->inode.max_limit = a->inode.max_limit; 2736 if (a->inode.size){ 2737 ierr = PetscMalloc((m+1)*sizeof(int),&c->inode.size);CHKERRQ(ierr); 2738 c->inode.node_count = a->inode.node_count; 2739 ierr = PetscMemcpy(c->inode.size,a->inode.size,(m+1)*sizeof(int));CHKERRQ(ierr); 2740 } else { 2741 c->inode.size = 0; 2742 c->inode.node_count = 0; 2743 } 2744 c->nz = a->nz; 2745 c->maxnz = a->maxnz; 2746 c->solve_work = 0; 2747 c->spptr = 0; /* Dangerous -I'm throwing away a->spptr */ 2748 C->preallocated = PETSC_TRUE; 2749 2750 *B = C; 2751 ierr = PetscFListDuplicate(A->qlist,&C->qlist);CHKERRQ(ierr); 2752 PetscFunctionReturn(0); 2753 } 2754 2755 EXTERN_C_BEGIN 2756 #undef __FUNCT__ 2757 #define __FUNCT__ "MatLoad_SeqAIJ" 2758 int MatLoad_SeqAIJ(PetscViewer viewer,MatType type,Mat *A) 2759 { 2760 Mat_SeqAIJ *a; 2761 Mat B; 2762 int i,nz,ierr,fd,header[4],size,*rowlengths = 0,M,N,shift; 2763 MPI_Comm comm; 2764 2765 PetscFunctionBegin; 2766 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 2767 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2768 if (size > 1) SETERRQ(PETSC_ERR_ARG_SIZ,"view must have one processor"); 2769 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 2770 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 2771 if (header[0] != MAT_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 2772 M = header[1]; N = header[2]; nz = header[3]; 2773 2774 if (nz < 0) { 2775 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 2776 } 2777 2778 /* read in row lengths */ 2779 ierr = PetscMalloc(M*sizeof(int),&rowlengths);CHKERRQ(ierr); 2780 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 2781 2782 /* create our matrix */ 2783 ierr = MatCreateSeqAIJ(comm,M,N,0,rowlengths,A);CHKERRQ(ierr); 2784 B = *A; 2785 a = (Mat_SeqAIJ*)B->data; 2786 shift = a->indexshift; 2787 2788 /* read in column indices and adjust for Fortran indexing*/ 2789 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 2790 if (shift) { 2791 for (i=0; i<nz; i++) { 2792 a->j[i] += 1; 2793 } 2794 } 2795 2796 /* read in nonzero values */ 2797 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 2798 2799 /* set matrix "i" values */ 2800 a->i[0] = -shift; 2801 for (i=1; i<= M; i++) { 2802 a->i[i] = a->i[i-1] + rowlengths[i-1]; 2803 a->ilen[i-1] = rowlengths[i-1]; 2804 } 2805 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 2806 2807 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2808 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2809 PetscFunctionReturn(0); 2810 } 2811 EXTERN_C_END 2812 2813 #undef __FUNCT__ 2814 #define __FUNCT__ "MatEqual_SeqAIJ" 2815 int MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg) 2816 { 2817 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data,*b = (Mat_SeqAIJ *)B->data; 2818 int ierr; 2819 PetscTruth flag; 2820 2821 PetscFunctionBegin; 2822 ierr = PetscTypeCompare((PetscObject)B,MATSEQAIJ,&flag);CHKERRQ(ierr); 2823 if (!flag) SETERRQ(PETSC_ERR_ARG_INCOMP,"Matrices must be same type"); 2824 2825 /* If the matrix dimensions are not equal,or no of nonzeros or shift */ 2826 if ((A->m != B->m) || (A->n != B->n) ||(a->nz != b->nz)|| (a->indexshift != b->indexshift)) { 2827 *flg = PETSC_FALSE; 2828 PetscFunctionReturn(0); 2829 } 2830 2831 /* if the a->i are the same */ 2832 ierr = PetscMemcmp(a->i,b->i,(A->m+1)*sizeof(int),flg);CHKERRQ(ierr); 2833 if (*flg == PETSC_FALSE) PetscFunctionReturn(0); 2834 2835 /* if a->j are the same */ 2836 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(int),flg);CHKERRQ(ierr); 2837 if (*flg == PETSC_FALSE) PetscFunctionReturn(0); 2838 2839 /* if a->a are the same */ 2840 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 2841 2842 PetscFunctionReturn(0); 2843 2844 } 2845 2846 #undef __FUNCT__ 2847 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 2848 /*@C 2849 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 2850 provided by the user. 2851 2852 Coolective on MPI_Comm 2853 2854 Input Parameters: 2855 + comm - must be an MPI communicator of size 1 2856 . m - number of rows 2857 . n - number of columns 2858 . i - row indices 2859 . j - column indices 2860 - a - matrix values 2861 2862 Output Parameter: 2863 . mat - the matrix 2864 2865 Level: intermediate 2866 2867 Notes: 2868 The i, j, and a arrays are not copied by this routine, the user must free these arrays 2869 once the matrix is destroyed 2870 2871 You cannot set new nonzero locations into this matrix, that will generate an error. 2872 2873 The i and j indices can be either 0- or 1 based 2874 2875 .seealso: MatCreate(), MatCreateMPIAIJ(), MatCreateSeqAIJ() 2876 2877 @*/ 2878 int MatCreateSeqAIJWithArrays(MPI_Comm comm,int m,int n,int* i,int*j,PetscScalar *a,Mat *mat) 2879 { 2880 int ierr,ii; 2881 Mat_SeqAIJ *aij; 2882 2883 PetscFunctionBegin; 2884 ierr = MatCreateSeqAIJ(comm,m,n,0,0,mat);CHKERRQ(ierr); 2885 aij = (Mat_SeqAIJ*)(*mat)->data; 2886 ierr = PetscFree(aij->a);CHKERRQ(ierr); 2887 2888 if (i[0] == 1) { 2889 aij->indexshift = -1; 2890 } else if (i[0]) { 2891 SETERRQ(1,"i (row indices) do not start with 0 or 1"); 2892 } 2893 aij->i = i; 2894 aij->j = j; 2895 aij->a = a; 2896 aij->singlemalloc = PETSC_FALSE; 2897 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 2898 aij->freedata = PETSC_FALSE; 2899 2900 for (ii=0; ii<m; ii++) { 2901 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 2902 #if defined(PETSC_USE_BOPT_g) 2903 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]); 2904 #endif 2905 } 2906 #if defined(PETSC_USE_BOPT_g) 2907 for (ii=0; ii<aij->i[m]; ii++) { 2908 if (j[ii] < -aij->indexshift) SETERRQ2(1,"Negative column index at location = %d index = %d",ii,j[ii]); 2909 if (j[ii] > n - 1 -aij->indexshift) SETERRQ2(1,"Column index to large at location = %d index = %d",ii,j[ii]); 2910 } 2911 #endif 2912 2913 /* changes indices to start at 0 */ 2914 if (i[0]) { 2915 aij->indexshift = 0; 2916 for (ii=0; ii<m; ii++) { 2917 i[ii]--; 2918 } 2919 for (ii=0; ii<i[m]; ii++) { 2920 j[ii]--; 2921 } 2922 } 2923 2924 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2925 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2926 PetscFunctionReturn(0); 2927 } 2928 2929 #undef __FUNCT__ 2930 #define __FUNCT__ "MatSetColoring_SeqAIJ" 2931 int MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 2932 { 2933 int ierr; 2934 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2935 2936 PetscFunctionBegin; 2937 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 2938 a->coloring = coloring; 2939 PetscFunctionReturn(0); 2940 } 2941 2942 #if defined(PETSC_HAVE_ADIC) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE) 2943 EXTERN_C_BEGIN 2944 #include "adic_utils.h" 2945 EXTERN_C_END 2946 2947 #undef __FUNCT__ 2948 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 2949 int MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 2950 { 2951 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2952 int m = A->m,*ii = a->i,*jj = a->j,nz,i,*color,j,nlen; 2953 PetscScalar *v = a->a,*values; 2954 char *cadvalues = (char *)advalues; 2955 2956 PetscFunctionBegin; 2957 if (!a->coloring) SETERRQ(1,"Coloring not set for matrix"); 2958 nlen = my_AD_GetDerivTypeSize(); 2959 color = a->coloring->colors; 2960 /* loop over rows */ 2961 for (i=0; i<m; i++) { 2962 nz = ii[i+1] - ii[i]; 2963 /* loop over columns putting computed value into matrix */ 2964 values = my_AD_GetGradArray(cadvalues); 2965 for (j=0; j<nz; j++) { 2966 *v++ = values[color[*jj++]]; 2967 } 2968 cadvalues += nlen; /* jump to next row of derivatives */ 2969 } 2970 PetscFunctionReturn(0); 2971 } 2972 2973 #else 2974 2975 #undef __FUNCT__ 2976 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 2977 int MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 2978 { 2979 PetscFunctionBegin; 2980 SETERRQ(1,"PETSc installed without ADIC"); 2981 } 2982 2983 #endif 2984 2985 #undef __FUNCT__ 2986 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 2987 int MatSetValuesAdifor_SeqAIJ(Mat A,int nl,void *advalues) 2988 { 2989 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2990 int m = A->m,*ii = a->i,*jj = a->j,nz,i,*color,j; 2991 PetscScalar *v = a->a,*values = (PetscScalar *)advalues; 2992 2993 PetscFunctionBegin; 2994 if (!a->coloring) SETERRQ(1,"Coloring not set for matrix"); 2995 color = a->coloring->colors; 2996 /* loop over rows */ 2997 for (i=0; i<m; i++) { 2998 nz = ii[i+1] - ii[i]; 2999 /* loop over columns putting computed value into matrix */ 3000 for (j=0; j<nz; j++) { 3001 *v++ = values[color[*jj++]]; 3002 } 3003 values += nl; /* jump to next row of derivatives */ 3004 } 3005 PetscFunctionReturn(0); 3006 } 3007 3008