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