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