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