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