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,sorted = a->sorted; 119 PetscInt *imax = a->imax,*ai = a->i,*ailen = a->ilen; 120 PetscErrorCode ierr; 121 PetscInt *aj = a->j,nonew = a->nonew; 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 (!sorted) low = 0; high = nrow; 150 while (high-low > 5) { 151 t = (low+high)/2; 152 if (rp[t] > col) high = t; 153 else low = t; 154 } 155 for (i=low; i<high; i++) { 156 if (rp[i] > col) break; 157 if (rp[i] == col) { 158 if (is == ADD_VALUES) ap[i] += value; 159 else ap[i] = value; 160 goto noinsert; 161 } 162 } 163 if (value == 0.0 && ignorezeroentries) goto noinsert; 164 if (nonew == 1) goto noinsert; 165 else if (nonew == -1) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at (%D,%D) in the matrix",row,col); 166 if (nrow >= rmax) { 167 /* there is no extra room in row, therefore enlarge */ 168 PetscInt new_nz = ai[A->m] + CHUNKSIZE,*new_i,*new_j; 169 size_t len; 170 PetscScalar *new_a; 171 172 if (nonew == -2) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at (%D,%D) in the matrix requiring new malloc()",row,col); 173 174 /* malloc new storage space */ 175 len = ((size_t) new_nz)*(sizeof(PetscInt)+sizeof(PetscScalar))+(A->m+1)*sizeof(PetscInt); 176 ierr = PetscMalloc(len,&new_a);CHKERRQ(ierr); 177 new_j = (PetscInt*)(new_a + new_nz); 178 new_i = new_j + new_nz; 179 180 /* copy over old data into new slots */ 181 for (ii=0; ii<row+1; ii++) {new_i[ii] = ai[ii];} 182 for (ii=row+1; ii<A->m+1; ii++) {new_i[ii] = ai[ii]+CHUNKSIZE;} 183 ierr = PetscMemcpy(new_j,aj,(ai[row]+nrow)*sizeof(PetscInt));CHKERRQ(ierr); 184 len = (((size_t) new_nz) - CHUNKSIZE - ai[row] - nrow ); 185 ierr = PetscMemcpy(new_j+ai[row]+nrow+CHUNKSIZE,aj+ai[row]+nrow,len*sizeof(PetscInt));CHKERRQ(ierr); 186 ierr = PetscMemcpy(new_a,aa,(((size_t) ai[row])+nrow)*sizeof(PetscScalar));CHKERRQ(ierr); 187 ierr = PetscMemcpy(new_a+ai[row]+nrow+CHUNKSIZE,aa+ai[row]+nrow,len*sizeof(PetscScalar));CHKERRQ(ierr); 188 /* free up old matrix storage */ 189 ierr = PetscFree(a->a);CHKERRQ(ierr); 190 if (!a->singlemalloc) { 191 ierr = PetscFree(a->i);CHKERRQ(ierr); 192 ierr = PetscFree(a->j);CHKERRQ(ierr); 193 } 194 aa = a->a = new_a; ai = a->i = new_i; aj = a->j = new_j; 195 a->singlemalloc = PETSC_TRUE; 196 197 rp = aj + ai[row]; ap = aa + ai[row] ; 198 rmax = imax[row] = imax[row] + CHUNKSIZE; 199 PetscLogObjectMemory(A,CHUNKSIZE*(sizeof(PetscInt) + sizeof(PetscScalar))); 200 a->maxnz += CHUNKSIZE; 201 a->reallocs++; 202 } 203 N = nrow++ - 1; a->nz++; 204 /* shift up all the later entries in this row */ 205 for (ii=N; ii>=i; ii--) { 206 rp[ii+1] = rp[ii]; 207 ap[ii+1] = ap[ii]; 208 } 209 rp[i] = col; 210 ap[i] = value; 211 noinsert:; 212 low = i + 1; 213 } 214 ailen[row] = nrow; 215 } 216 A->same_nonzero = PETSC_FALSE; 217 PetscFunctionReturn(0); 218 } 219 220 #undef __FUNCT__ 221 #define __FUNCT__ "MatGetValues_SeqAIJ" 222 PetscErrorCode MatGetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],PetscScalar v[]) 223 { 224 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 225 PetscInt *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j; 226 PetscInt *ai = a->i,*ailen = a->ilen; 227 PetscScalar *ap,*aa = a->a,zero = 0.0; 228 229 PetscFunctionBegin; 230 for (k=0; k<m; k++) { /* loop over rows */ 231 row = im[k]; 232 if (row < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Negative row: %D",row); 233 if (row >= A->m) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->m-1); 234 rp = aj + ai[row]; ap = aa + ai[row]; 235 nrow = ailen[row]; 236 for (l=0; l<n; l++) { /* loop over columns */ 237 if (in[l] < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Negative column: %D",in[l]); 238 if (in[l] >= A->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->n-1); 239 col = in[l] ; 240 high = nrow; low = 0; /* assume unsorted */ 241 while (high-low > 5) { 242 t = (low+high)/2; 243 if (rp[t] > col) high = t; 244 else low = t; 245 } 246 for (i=low; i<high; i++) { 247 if (rp[i] > col) break; 248 if (rp[i] == col) { 249 *v++ = ap[i]; 250 goto finished; 251 } 252 } 253 *v++ = zero; 254 finished:; 255 } 256 } 257 PetscFunctionReturn(0); 258 } 259 260 261 #undef __FUNCT__ 262 #define __FUNCT__ "MatView_SeqAIJ_Binary" 263 PetscErrorCode MatView_SeqAIJ_Binary(Mat A,PetscViewer viewer) 264 { 265 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 266 PetscErrorCode ierr; 267 PetscInt i,*col_lens; 268 int fd; 269 270 PetscFunctionBegin; 271 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 272 ierr = PetscMalloc((4+A->m)*sizeof(PetscInt),&col_lens);CHKERRQ(ierr); 273 col_lens[0] = MAT_FILE_COOKIE; 274 col_lens[1] = A->m; 275 col_lens[2] = A->n; 276 col_lens[3] = a->nz; 277 278 /* store lengths of each row and write (including header) to file */ 279 for (i=0; i<A->m; i++) { 280 col_lens[4+i] = a->i[i+1] - a->i[i]; 281 } 282 ierr = PetscBinaryWrite(fd,col_lens,4+A->m,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 283 ierr = PetscFree(col_lens);CHKERRQ(ierr); 284 285 /* store column indices (zero start index) */ 286 ierr = PetscBinaryWrite(fd,a->j,a->nz,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 287 288 /* store nonzero values */ 289 ierr = PetscBinaryWrite(fd,a->a,a->nz,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 290 PetscFunctionReturn(0); 291 } 292 293 EXTERN PetscErrorCode MatSeqAIJFactorInfo_Matlab(Mat,PetscViewer); 294 295 #undef __FUNCT__ 296 #define __FUNCT__ "MatView_SeqAIJ_ASCII" 297 PetscErrorCode MatView_SeqAIJ_ASCII(Mat A,PetscViewer viewer) 298 { 299 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 300 PetscErrorCode ierr; 301 PetscInt i,j,m = A->m,shift=0; 302 char *name; 303 PetscViewerFormat format; 304 305 PetscFunctionBegin; 306 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 307 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 308 if (format == PETSC_VIEWER_ASCII_INFO_DETAIL || format == PETSC_VIEWER_ASCII_INFO) { 309 if (a->inode.size) { 310 ierr = PetscViewerASCIIPrintf(viewer,"using I-node routines: found %D nodes, limit used is %D\n",a->inode.node_count,a->inode.limit);CHKERRQ(ierr); 311 } else { 312 ierr = PetscViewerASCIIPrintf(viewer,"not using I-node routines\n");CHKERRQ(ierr); 313 } 314 } else if (format == PETSC_VIEWER_ASCII_MATLAB) { 315 PetscInt nofinalvalue = 0; 316 if ((a->i[m] == a->i[m-1]) || (a->j[a->nz-1] != A->n-!shift)) { 317 nofinalvalue = 1; 318 } 319 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 320 ierr = PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",m,A->n);CHKERRQ(ierr); 321 ierr = PetscViewerASCIIPrintf(viewer,"%% Nonzeros = %D \n",a->nz);CHKERRQ(ierr); 322 ierr = PetscViewerASCIIPrintf(viewer,"zzz = zeros(%D,3);\n",a->nz+nofinalvalue);CHKERRQ(ierr); 323 ierr = PetscViewerASCIIPrintf(viewer,"zzz = [\n");CHKERRQ(ierr); 324 325 for (i=0; i<m; i++) { 326 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 327 #if defined(PETSC_USE_COMPLEX) 328 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); 329 #else 330 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %18.16e\n",i+1,a->j[j]+!shift,a->a[j]);CHKERRQ(ierr); 331 #endif 332 } 333 } 334 if (nofinalvalue) { 335 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %18.16e\n",m,A->n,0.0);CHKERRQ(ierr); 336 } 337 ierr = PetscViewerASCIIPrintf(viewer,"];\n %s = spconvert(zzz);\n",name);CHKERRQ(ierr); 338 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 339 } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) { 340 PetscFunctionReturn(0); 341 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 342 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 343 for (i=0; i<m; i++) { 344 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 345 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 346 #if defined(PETSC_USE_COMPLEX) 347 if (PetscImaginaryPart(a->a[j]) > 0.0 && PetscRealPart(a->a[j]) != 0.0) { 348 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 349 } else if (PetscImaginaryPart(a->a[j]) < 0.0 && PetscRealPart(a->a[j]) != 0.0) { 350 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g - %g i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 351 } else if (PetscRealPart(a->a[j]) != 0.0) { 352 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 353 } 354 #else 355 if (a->a[j] != 0.0) {ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr);} 356 #endif 357 } 358 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 359 } 360 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 361 } else if (format == PETSC_VIEWER_ASCII_SYMMODU) { 362 PetscInt nzd=0,fshift=1,*sptr; 363 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 364 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&sptr);CHKERRQ(ierr); 365 for (i=0; i<m; i++) { 366 sptr[i] = nzd+1; 367 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 368 if (a->j[j] >= i) { 369 #if defined(PETSC_USE_COMPLEX) 370 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) nzd++; 371 #else 372 if (a->a[j] != 0.0) nzd++; 373 #endif 374 } 375 } 376 } 377 sptr[m] = nzd+1; 378 ierr = PetscViewerASCIIPrintf(viewer," %D %D\n\n",m,nzd);CHKERRQ(ierr); 379 for (i=0; i<m+1; i+=6) { 380 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);} 381 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);} 382 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);} 383 else if (i+1<m) {ierr = PetscViewerASCIIPrintf(viewer," %D %D %D\n",sptr[i],sptr[i+1],sptr[i+2]);CHKERRQ(ierr);} 384 else if (i<m) {ierr = PetscViewerASCIIPrintf(viewer," %D %D\n",sptr[i],sptr[i+1]);CHKERRQ(ierr);} 385 else {ierr = PetscViewerASCIIPrintf(viewer," %D\n",sptr[i]);CHKERRQ(ierr);} 386 } 387 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 388 ierr = PetscFree(sptr);CHKERRQ(ierr); 389 for (i=0; i<m; i++) { 390 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 391 if (a->j[j] >= i) {ierr = PetscViewerASCIIPrintf(viewer," %D ",a->j[j]+fshift);CHKERRQ(ierr);} 392 } 393 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 394 } 395 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 396 for (i=0; i<m; i++) { 397 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 398 if (a->j[j] >= i) { 399 #if defined(PETSC_USE_COMPLEX) 400 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) { 401 ierr = PetscViewerASCIIPrintf(viewer," %18.16e %18.16e ",PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 402 } 403 #else 404 if (a->a[j] != 0.0) {ierr = PetscViewerASCIIPrintf(viewer," %18.16e ",a->a[j]);CHKERRQ(ierr);} 405 #endif 406 } 407 } 408 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 409 } 410 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 411 } else if (format == PETSC_VIEWER_ASCII_DENSE) { 412 PetscInt cnt = 0,jcnt; 413 PetscScalar value; 414 415 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 416 for (i=0; i<m; i++) { 417 jcnt = 0; 418 for (j=0; j<A->n; j++) { 419 if (jcnt < a->i[i+1]-a->i[i] && j == a->j[cnt]) { 420 value = a->a[cnt++]; 421 jcnt++; 422 } else { 423 value = 0.0; 424 } 425 #if defined(PETSC_USE_COMPLEX) 426 ierr = PetscViewerASCIIPrintf(viewer," %7.5e+%7.5e i ",PetscRealPart(value),PetscImaginaryPart(value));CHKERRQ(ierr); 427 #else 428 ierr = PetscViewerASCIIPrintf(viewer," %7.5e ",value);CHKERRQ(ierr); 429 #endif 430 } 431 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 432 } 433 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 434 } else { 435 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 436 for (i=0; i<m; i++) { 437 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 438 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 439 #if defined(PETSC_USE_COMPLEX) 440 if (PetscImaginaryPart(a->a[j]) > 0.0) { 441 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g + %g i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 442 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 443 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g - %g i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 444 } else { 445 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 446 } 447 #else 448 ierr = PetscViewerASCIIPrintf(viewer," (%D, %g) ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr); 449 #endif 450 } 451 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 452 } 453 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 454 } 455 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 456 PetscFunctionReturn(0); 457 } 458 459 #undef __FUNCT__ 460 #define __FUNCT__ "MatView_SeqAIJ_Draw_Zoom" 461 PetscErrorCode MatView_SeqAIJ_Draw_Zoom(PetscDraw draw,void *Aa) 462 { 463 Mat A = (Mat) Aa; 464 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 465 PetscErrorCode ierr; 466 PetscInt i,j,m = A->m,color; 467 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r,maxv = 0.0; 468 PetscViewer viewer; 469 PetscViewerFormat format; 470 471 PetscFunctionBegin; 472 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 473 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 474 475 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 476 /* loop over matrix elements drawing boxes */ 477 478 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 479 /* Blue for negative, Cyan for zero and Red for positive */ 480 color = PETSC_DRAW_BLUE; 481 for (i=0; i<m; i++) { 482 y_l = m - i - 1.0; y_r = y_l + 1.0; 483 for (j=a->i[i]; j<a->i[i+1]; j++) { 484 x_l = a->j[j] ; x_r = x_l + 1.0; 485 #if defined(PETSC_USE_COMPLEX) 486 if (PetscRealPart(a->a[j]) >= 0.) continue; 487 #else 488 if (a->a[j] >= 0.) continue; 489 #endif 490 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 491 } 492 } 493 color = PETSC_DRAW_CYAN; 494 for (i=0; i<m; i++) { 495 y_l = m - i - 1.0; y_r = y_l + 1.0; 496 for (j=a->i[i]; j<a->i[i+1]; j++) { 497 x_l = a->j[j]; x_r = x_l + 1.0; 498 if (a->a[j] != 0.) continue; 499 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 500 } 501 } 502 color = PETSC_DRAW_RED; 503 for (i=0; i<m; i++) { 504 y_l = m - i - 1.0; y_r = y_l + 1.0; 505 for (j=a->i[i]; j<a->i[i+1]; j++) { 506 x_l = a->j[j]; x_r = x_l + 1.0; 507 #if defined(PETSC_USE_COMPLEX) 508 if (PetscRealPart(a->a[j]) <= 0.) continue; 509 #else 510 if (a->a[j] <= 0.) continue; 511 #endif 512 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 513 } 514 } 515 } else { 516 /* use contour shading to indicate magnitude of values */ 517 /* first determine max of all nonzero values */ 518 PetscInt nz = a->nz,count; 519 PetscDraw popup; 520 PetscReal scale; 521 522 for (i=0; i<nz; i++) { 523 if (PetscAbsScalar(a->a[i]) > maxv) maxv = PetscAbsScalar(a->a[i]); 524 } 525 scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/maxv; 526 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 527 if (popup) {ierr = PetscDrawScalePopup(popup,0.0,maxv);CHKERRQ(ierr);} 528 count = 0; 529 for (i=0; i<m; i++) { 530 y_l = m - i - 1.0; y_r = y_l + 1.0; 531 for (j=a->i[i]; j<a->i[i+1]; j++) { 532 x_l = a->j[j]; x_r = x_l + 1.0; 533 color = PETSC_DRAW_BASIC_COLORS + (PetscInt)(scale*PetscAbsScalar(a->a[count])); 534 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 535 count++; 536 } 537 } 538 } 539 PetscFunctionReturn(0); 540 } 541 542 #undef __FUNCT__ 543 #define __FUNCT__ "MatView_SeqAIJ_Draw" 544 PetscErrorCode MatView_SeqAIJ_Draw(Mat A,PetscViewer viewer) 545 { 546 PetscErrorCode ierr; 547 PetscDraw draw; 548 PetscReal xr,yr,xl,yl,h,w; 549 PetscTruth isnull; 550 551 PetscFunctionBegin; 552 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 553 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 554 if (isnull) PetscFunctionReturn(0); 555 556 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 557 xr = A->n; yr = A->m; h = yr/10.0; w = xr/10.0; 558 xr += w; yr += h; xl = -w; yl = -h; 559 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 560 ierr = PetscDrawZoom(draw,MatView_SeqAIJ_Draw_Zoom,A);CHKERRQ(ierr); 561 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);CHKERRQ(ierr); 562 PetscFunctionReturn(0); 563 } 564 565 #undef __FUNCT__ 566 #define __FUNCT__ "MatView_SeqAIJ" 567 PetscErrorCode MatView_SeqAIJ(Mat A,PetscViewer viewer) 568 { 569 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 570 PetscErrorCode ierr; 571 PetscTruth issocket,iascii,isbinary,isdraw; 572 573 PetscFunctionBegin; 574 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_SOCKET,&issocket);CHKERRQ(ierr); 575 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 576 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_BINARY,&isbinary);CHKERRQ(ierr); 577 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_DRAW,&isdraw);CHKERRQ(ierr); 578 if (issocket) { 579 ierr = PetscViewerSocketPutSparse_Private(viewer,A->m,A->n,a->nz,a->a,a->i,a->j);CHKERRQ(ierr); 580 } else if (iascii) { 581 ierr = MatView_SeqAIJ_ASCII(A,viewer);CHKERRQ(ierr); 582 } else if (isbinary) { 583 ierr = MatView_SeqAIJ_Binary(A,viewer);CHKERRQ(ierr); 584 } else if (isdraw) { 585 ierr = MatView_SeqAIJ_Draw(A,viewer);CHKERRQ(ierr); 586 } else { 587 SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported by SeqAIJ matrices",((PetscObject)viewer)->type_name); 588 } 589 PetscFunctionReturn(0); 590 } 591 592 #undef __FUNCT__ 593 #define __FUNCT__ "MatAssemblyEnd_SeqAIJ" 594 PetscErrorCode MatAssemblyEnd_SeqAIJ(Mat A,MatAssemblyType mode) 595 { 596 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 597 PetscErrorCode ierr; 598 PetscInt fshift = 0,i,j,*ai = a->i,*aj = a->j,*imax = a->imax; 599 PetscInt m = A->m,*ip,N,*ailen = a->ilen,rmax = 0; 600 PetscScalar *aa = a->a,*ap; 601 PetscReal ratio=0.6; 602 603 PetscFunctionBegin; 604 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 605 606 if (m) rmax = ailen[0]; /* determine row with most nonzeros */ 607 for (i=1; i<m; i++) { 608 /* move each row back by the amount of empty slots (fshift) before it*/ 609 fshift += imax[i-1] - ailen[i-1]; 610 rmax = PetscMax(rmax,ailen[i]); 611 if (fshift) { 612 ip = aj + ai[i] ; 613 ap = aa + ai[i] ; 614 N = ailen[i]; 615 for (j=0; j<N; j++) { 616 ip[j-fshift] = ip[j]; 617 ap[j-fshift] = ap[j]; 618 } 619 } 620 ai[i] = ai[i-1] + ailen[i-1]; 621 } 622 if (m) { 623 fshift += imax[m-1] - ailen[m-1]; 624 ai[m] = ai[m-1] + ailen[m-1]; 625 } 626 /* reset ilen and imax for each row */ 627 for (i=0; i<m; i++) { 628 ailen[i] = imax[i] = ai[i+1] - ai[i]; 629 } 630 a->nz = ai[m]; 631 632 /* diagonals may have moved, so kill the diagonal pointers */ 633 if (fshift && a->diag) { 634 ierr = PetscFree(a->diag);CHKERRQ(ierr); 635 PetscLogObjectMemory(A,-(m+1)*sizeof(PetscInt)); 636 a->diag = 0; 637 } 638 PetscLogInfo(A,"MatAssemblyEnd_SeqAIJ:Matrix size: %D X %D; storage space: %D unneeded,%D used\n",m,A->n,fshift,a->nz); 639 PetscLogInfo(A,"MatAssemblyEnd_SeqAIJ:Number of mallocs during MatSetValues() is %D\n",a->reallocs); 640 PetscLogInfo(A,"MatAssemblyEnd_SeqAIJ:Maximum nonzeros in any row is %D\n",rmax); 641 a->reallocs = 0; 642 A->info.nz_unneeded = (double)fshift; 643 a->rmax = rmax; 644 645 /* check for identical nodes. If found, use inode functions */ 646 if (!a->inode.checked){ 647 ierr = Mat_AIJ_CheckInode(A,(PetscTruth)(!fshift));CHKERRQ(ierr); 648 } 649 650 /* check for zero rows. If found a large number of zero rows, use CompressedRow functions */ 651 if (!a->inode.use && a->compressedrow.use && !A->same_nonzero){ 652 ierr = Mat_CheckCompressedRow(A,&a->compressedrow,a->i,ratio);CHKERRQ(ierr); 653 } 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 PetscLogObjectMemory(A,(m+1)*sizeof(PetscInt)); 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 PetscLogObjectParent(inA,a->icol); 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,&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 PetscLogObjectParent(coloring,coloring->w1); 1998 ierr = VecDuplicate(x1,&coloring->w2);CHKERRQ(ierr); 1999 PetscLogObjectParent(coloring,coloring->w2); 2000 ierr = VecDuplicate(x1,&coloring->w3);CHKERRQ(ierr); 2001 PetscLogObjectParent(coloring,coloring->w3); 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_SeqAIJ_SeqAIJ, 2272 MatPtAPSymbolic_SeqAIJ_SeqAIJ, 2273 /*95*/ MatPtAPNumeric_SeqAIJ_SeqAIJ, 2274 MatMatMultTranspose_SeqAIJ_SeqAIJ, 2275 MatMatMultTransposeSymbolic_SeqAIJ_SeqAIJ, 2276 MatMatMultTransposeNumeric_SeqAIJ_SeqAIJ, 2277 }; 2278 2279 EXTERN_C_BEGIN 2280 #undef __FUNCT__ 2281 #define __FUNCT__ "MatSeqAIJSetColumnIndices_SeqAIJ" 2282 PetscErrorCode MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 2283 { 2284 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2285 PetscInt i,nz,n; 2286 2287 PetscFunctionBegin; 2288 2289 nz = aij->maxnz; 2290 n = mat->n; 2291 for (i=0; i<nz; i++) { 2292 aij->j[i] = indices[i]; 2293 } 2294 aij->nz = nz; 2295 for (i=0; i<n; i++) { 2296 aij->ilen[i] = aij->imax[i]; 2297 } 2298 2299 PetscFunctionReturn(0); 2300 } 2301 EXTERN_C_END 2302 2303 #undef __FUNCT__ 2304 #define __FUNCT__ "MatSeqAIJSetColumnIndices" 2305 /*@ 2306 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 2307 in the matrix. 2308 2309 Input Parameters: 2310 + mat - the SeqAIJ matrix 2311 - indices - the column indices 2312 2313 Level: advanced 2314 2315 Notes: 2316 This can be called if you have precomputed the nonzero structure of the 2317 matrix and want to provide it to the matrix object to improve the performance 2318 of the MatSetValues() operation. 2319 2320 You MUST have set the correct numbers of nonzeros per row in the call to 2321 MatCreateSeqAIJ(). 2322 2323 MUST be called before any calls to MatSetValues(); 2324 2325 The indices should start with zero, not one. 2326 2327 @*/ 2328 PetscErrorCode MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 2329 { 2330 PetscErrorCode ierr,(*f)(Mat,PetscInt *); 2331 2332 PetscFunctionBegin; 2333 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 2334 PetscValidPointer(indices,2); 2335 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqAIJSetColumnIndices_C",(void (**)(void))&f);CHKERRQ(ierr); 2336 if (f) { 2337 ierr = (*f)(mat,indices);CHKERRQ(ierr); 2338 } else { 2339 SETERRQ(PETSC_ERR_SUP,"Wrong type of matrix to set column indices"); 2340 } 2341 PetscFunctionReturn(0); 2342 } 2343 2344 /* ----------------------------------------------------------------------------------------*/ 2345 2346 EXTERN_C_BEGIN 2347 #undef __FUNCT__ 2348 #define __FUNCT__ "MatStoreValues_SeqAIJ" 2349 PetscErrorCode MatStoreValues_SeqAIJ(Mat mat) 2350 { 2351 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2352 PetscErrorCode ierr; 2353 size_t nz = aij->i[mat->m]; 2354 2355 PetscFunctionBegin; 2356 if (aij->nonew != 1) { 2357 SETERRQ(PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 2358 } 2359 2360 /* allocate space for values if not already there */ 2361 if (!aij->saved_values) { 2362 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 2363 } 2364 2365 /* copy values over */ 2366 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2367 PetscFunctionReturn(0); 2368 } 2369 EXTERN_C_END 2370 2371 #undef __FUNCT__ 2372 #define __FUNCT__ "MatStoreValues" 2373 /*@ 2374 MatStoreValues - Stashes a copy of the matrix values; this allows, for 2375 example, reuse of the linear part of a Jacobian, while recomputing the 2376 nonlinear portion. 2377 2378 Collect on Mat 2379 2380 Input Parameters: 2381 . mat - the matrix (currently on AIJ matrices support this option) 2382 2383 Level: advanced 2384 2385 Common Usage, with SNESSolve(): 2386 $ Create Jacobian matrix 2387 $ Set linear terms into matrix 2388 $ Apply boundary conditions to matrix, at this time matrix must have 2389 $ final nonzero structure (i.e. setting the nonlinear terms and applying 2390 $ boundary conditions again will not change the nonzero structure 2391 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 2392 $ ierr = MatStoreValues(mat); 2393 $ Call SNESSetJacobian() with matrix 2394 $ In your Jacobian routine 2395 $ ierr = MatRetrieveValues(mat); 2396 $ Set nonlinear terms in matrix 2397 2398 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 2399 $ // build linear portion of Jacobian 2400 $ ierr = MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); 2401 $ ierr = MatStoreValues(mat); 2402 $ loop over nonlinear iterations 2403 $ ierr = MatRetrieveValues(mat); 2404 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 2405 $ // call MatAssemblyBegin/End() on matrix 2406 $ Solve linear system with Jacobian 2407 $ endloop 2408 2409 Notes: 2410 Matrix must already be assemblied before calling this routine 2411 Must set the matrix option MatSetOption(mat,MAT_NO_NEW_NONZERO_LOCATIONS); before 2412 calling this routine. 2413 2414 .seealso: MatRetrieveValues() 2415 2416 @*/ 2417 PetscErrorCode MatStoreValues(Mat mat) 2418 { 2419 PetscErrorCode ierr,(*f)(Mat); 2420 2421 PetscFunctionBegin; 2422 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 2423 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2424 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2425 2426 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatStoreValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2427 if (f) { 2428 ierr = (*f)(mat);CHKERRQ(ierr); 2429 } else { 2430 SETERRQ(PETSC_ERR_SUP,"Wrong type of matrix to store values"); 2431 } 2432 PetscFunctionReturn(0); 2433 } 2434 2435 EXTERN_C_BEGIN 2436 #undef __FUNCT__ 2437 #define __FUNCT__ "MatRetrieveValues_SeqAIJ" 2438 PetscErrorCode MatRetrieveValues_SeqAIJ(Mat mat) 2439 { 2440 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2441 PetscErrorCode ierr; 2442 PetscInt nz = aij->i[mat->m]; 2443 2444 PetscFunctionBegin; 2445 if (aij->nonew != 1) { 2446 SETERRQ(PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NO_NEW_NONZERO_LOCATIONS);first"); 2447 } 2448 if (!aij->saved_values) { 2449 SETERRQ(PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 2450 } 2451 /* copy values over */ 2452 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2453 PetscFunctionReturn(0); 2454 } 2455 EXTERN_C_END 2456 2457 #undef __FUNCT__ 2458 #define __FUNCT__ "MatRetrieveValues" 2459 /*@ 2460 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 2461 example, reuse of the linear part of a Jacobian, while recomputing the 2462 nonlinear portion. 2463 2464 Collect on Mat 2465 2466 Input Parameters: 2467 . mat - the matrix (currently on AIJ matrices support this option) 2468 2469 Level: advanced 2470 2471 .seealso: MatStoreValues() 2472 2473 @*/ 2474 PetscErrorCode MatRetrieveValues(Mat mat) 2475 { 2476 PetscErrorCode ierr,(*f)(Mat); 2477 2478 PetscFunctionBegin; 2479 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 2480 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2481 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2482 2483 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatRetrieveValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2484 if (f) { 2485 ierr = (*f)(mat);CHKERRQ(ierr); 2486 } else { 2487 SETERRQ(PETSC_ERR_SUP,"Wrong type of matrix to retrieve values"); 2488 } 2489 PetscFunctionReturn(0); 2490 } 2491 2492 2493 /* --------------------------------------------------------------------------------*/ 2494 #undef __FUNCT__ 2495 #define __FUNCT__ "MatCreateSeqAIJ" 2496 /*@C 2497 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 2498 (the default parallel PETSc format). For good matrix assembly performance 2499 the user should preallocate the matrix storage by setting the parameter nz 2500 (or the array nnz). By setting these parameters accurately, performance 2501 during matrix assembly can be increased by more than a factor of 50. 2502 2503 Collective on MPI_Comm 2504 2505 Input Parameters: 2506 + comm - MPI communicator, set to PETSC_COMM_SELF 2507 . m - number of rows 2508 . n - number of columns 2509 . nz - number of nonzeros per row (same for all rows) 2510 - nnz - array containing the number of nonzeros in the various rows 2511 (possibly different for each row) or PETSC_NULL 2512 2513 Output Parameter: 2514 . A - the matrix 2515 2516 Notes: 2517 If nnz is given then nz is ignored 2518 2519 The AIJ format (also called the Yale sparse matrix format or 2520 compressed row storage), is fully compatible with standard Fortran 77 2521 storage. That is, the stored row and column indices can begin at 2522 either one (as in Fortran) or zero. See the users' manual for details. 2523 2524 Specify the preallocated storage with either nz or nnz (not both). 2525 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2526 allocation. For large problems you MUST preallocate memory or you 2527 will get TERRIBLE performance, see the users' manual chapter on matrices. 2528 2529 By default, this format uses inodes (identical nodes) when possible, to 2530 improve numerical efficiency of matrix-vector products and solves. We 2531 search for consecutive rows with the same nonzero structure, thereby 2532 reusing matrix information to achieve increased efficiency. 2533 2534 Options Database Keys: 2535 + -mat_aij_no_inode - Do not use inodes 2536 . -mat_aij_inode_limit <limit> - Sets inode limit (max limit=5) 2537 - -mat_aij_oneindex - Internally use indexing starting at 1 2538 rather than 0. Note that when calling MatSetValues(), 2539 the user still MUST index entries starting at 0! 2540 2541 Level: intermediate 2542 2543 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2544 2545 @*/ 2546 PetscErrorCode MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 2547 { 2548 PetscErrorCode ierr; 2549 2550 PetscFunctionBegin; 2551 ierr = MatCreate(comm,m,n,m,n,A);CHKERRQ(ierr); 2552 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 2553 ierr = MatSeqAIJSetPreallocation(*A,nz,nnz);CHKERRQ(ierr); 2554 PetscFunctionReturn(0); 2555 } 2556 2557 #define SKIP_ALLOCATION -4 2558 2559 #undef __FUNCT__ 2560 #define __FUNCT__ "MatSeqAIJSetPreallocation" 2561 /*@C 2562 MatSeqAIJSetPreallocation - For good matrix assembly performance 2563 the user should preallocate the matrix storage by setting the parameter nz 2564 (or the array nnz). By setting these parameters accurately, performance 2565 during matrix assembly can be increased by more than a factor of 50. 2566 2567 Collective on MPI_Comm 2568 2569 Input Parameters: 2570 + comm - MPI communicator, set to PETSC_COMM_SELF 2571 . m - number of rows 2572 . n - number of columns 2573 . nz - number of nonzeros per row (same for all rows) 2574 - nnz - array containing the number of nonzeros in the various rows 2575 (possibly different for each row) or PETSC_NULL 2576 2577 Output Parameter: 2578 . A - the matrix 2579 2580 Notes: 2581 If nnz is given then nz is ignored 2582 2583 The AIJ format (also called the Yale sparse matrix format or 2584 compressed row storage), is fully compatible with standard Fortran 77 2585 storage. That is, the stored row and column indices can begin at 2586 either one (as in Fortran) or zero. See the users' manual for details. 2587 2588 Specify the preallocated storage with either nz or nnz (not both). 2589 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2590 allocation. For large problems you MUST preallocate memory or you 2591 will get TERRIBLE performance, see the users' manual chapter on matrices. 2592 2593 By default, this format uses inodes (identical nodes) when possible, to 2594 improve numerical efficiency of matrix-vector products and solves. We 2595 search for consecutive rows with the same nonzero structure, thereby 2596 reusing matrix information to achieve increased efficiency. 2597 2598 Options Database Keys: 2599 + -mat_aij_no_inode - Do not use inodes 2600 . -mat_aij_inode_limit <limit> - Sets inode limit (max limit=5) 2601 - -mat_aij_oneindex - Internally use indexing starting at 1 2602 rather than 0. Note that when calling MatSetValues(), 2603 the user still MUST index entries starting at 0! 2604 2605 Level: intermediate 2606 2607 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2608 2609 @*/ 2610 PetscErrorCode MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 2611 { 2612 PetscErrorCode ierr,(*f)(Mat,PetscInt,const PetscInt[]); 2613 2614 PetscFunctionBegin; 2615 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 2616 if (f) { 2617 ierr = (*f)(B,nz,nnz);CHKERRQ(ierr); 2618 } 2619 PetscFunctionReturn(0); 2620 } 2621 2622 EXTERN_C_BEGIN 2623 #undef __FUNCT__ 2624 #define __FUNCT__ "MatSeqAIJSetPreallocation_SeqAIJ" 2625 PetscErrorCode MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,PetscInt *nnz) 2626 { 2627 Mat_SeqAIJ *b; 2628 size_t len = 0; 2629 PetscTruth skipallocation = PETSC_FALSE; 2630 PetscErrorCode ierr; 2631 PetscInt i; 2632 2633 PetscFunctionBegin; 2634 2635 if (nz == SKIP_ALLOCATION) { 2636 skipallocation = PETSC_TRUE; 2637 nz = 0; 2638 } 2639 2640 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 2641 if (nz < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 2642 if (nnz) { 2643 for (i=0; i<B->m; i++) { 2644 if (nnz[i] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]); 2645 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); 2646 } 2647 } 2648 2649 B->preallocated = PETSC_TRUE; 2650 b = (Mat_SeqAIJ*)B->data; 2651 2652 ierr = PetscMalloc((B->m+1)*sizeof(PetscInt),&b->imax);CHKERRQ(ierr); 2653 if (!nnz) { 2654 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 2655 else if (nz <= 0) nz = 1; 2656 for (i=0; i<B->m; i++) b->imax[i] = nz; 2657 nz = nz*B->m; 2658 } else { 2659 nz = 0; 2660 for (i=0; i<B->m; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 2661 } 2662 2663 if (!skipallocation) { 2664 /* allocate the matrix space */ 2665 len = ((size_t) nz)*(sizeof(PetscInt) + sizeof(PetscScalar)) + (B->m+1)*sizeof(PetscInt); 2666 ierr = PetscMalloc(len,&b->a);CHKERRQ(ierr); 2667 b->j = (PetscInt*)(b->a + nz); 2668 ierr = PetscMemzero(b->j,nz*sizeof(PetscInt));CHKERRQ(ierr); 2669 b->i = b->j + nz; 2670 b->i[0] = 0; 2671 for (i=1; i<B->m+1; i++) { 2672 b->i[i] = b->i[i-1] + b->imax[i-1]; 2673 } 2674 b->singlemalloc = PETSC_TRUE; 2675 b->freedata = PETSC_TRUE; 2676 } else { 2677 b->freedata = PETSC_FALSE; 2678 } 2679 2680 /* b->ilen will count nonzeros in each row so far. */ 2681 ierr = PetscMalloc((B->m+1)*sizeof(PetscInt),&b->ilen);CHKERRQ(ierr); 2682 PetscLogObjectMemory(B,len+2*(B->m+1)*sizeof(PetscInt)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2683 for (i=0; i<B->m; i++) { b->ilen[i] = 0;} 2684 2685 b->nz = 0; 2686 b->maxnz = nz; 2687 B->info.nz_unneeded = (double)b->maxnz; 2688 PetscFunctionReturn(0); 2689 } 2690 EXTERN_C_END 2691 2692 /*MC 2693 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 2694 based on compressed sparse row format. 2695 2696 Options Database Keys: 2697 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 2698 2699 Level: beginner 2700 2701 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 2702 M*/ 2703 2704 EXTERN_C_BEGIN 2705 #undef __FUNCT__ 2706 #define __FUNCT__ "MatCreate_SeqAIJ" 2707 PetscErrorCode MatCreate_SeqAIJ(Mat B) 2708 { 2709 Mat_SeqAIJ *b; 2710 PetscErrorCode ierr; 2711 PetscMPIInt size; 2712 2713 PetscFunctionBegin; 2714 ierr = MPI_Comm_size(B->comm,&size);CHKERRQ(ierr); 2715 if (size > 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 2716 2717 B->m = B->M = PetscMax(B->m,B->M); 2718 B->n = B->N = PetscMax(B->n,B->N); 2719 2720 ierr = PetscNew(Mat_SeqAIJ,&b);CHKERRQ(ierr); 2721 B->data = (void*)b; 2722 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 2723 B->factor = 0; 2724 B->lupivotthreshold = 1.0; 2725 B->mapping = 0; 2726 ierr = PetscOptionsGetReal(B->prefix,"-mat_lu_pivotthreshold",&B->lupivotthreshold,PETSC_NULL);CHKERRQ(ierr); 2727 ierr = PetscOptionsHasName(B->prefix,"-pc_ilu_preserve_row_sums",&b->ilu_preserve_row_sums);CHKERRQ(ierr); 2728 b->row = 0; 2729 b->col = 0; 2730 b->icol = 0; 2731 b->reallocs = 0; 2732 b->lu_shift = PETSC_FALSE; 2733 2734 ierr = PetscMapCreateMPI(B->comm,B->m,B->m,&B->rmap);CHKERRQ(ierr); 2735 ierr = PetscMapCreateMPI(B->comm,B->n,B->n,&B->cmap);CHKERRQ(ierr); 2736 2737 b->sorted = PETSC_FALSE; 2738 b->ignorezeroentries = PETSC_FALSE; 2739 b->roworiented = PETSC_TRUE; 2740 b->nonew = 0; 2741 b->diag = 0; 2742 b->solve_work = 0; 2743 B->spptr = 0; 2744 b->inode.use = PETSC_TRUE; 2745 b->inode.node_count = 0; 2746 b->inode.size = 0; 2747 b->inode.limit = 5; 2748 b->inode.max_limit = 5; 2749 b->saved_values = 0; 2750 b->idiag = 0; 2751 b->ssor = 0; 2752 b->keepzeroedrows = PETSC_FALSE; 2753 b->xtoy = 0; 2754 b->XtoY = 0; 2755 b->compressedrow.use = PETSC_FALSE; 2756 b->compressedrow.nrows = B->m; 2757 b->compressedrow.i = PETSC_NULL; 2758 b->compressedrow.rindex = PETSC_NULL; 2759 b->compressedrow.checked = PETSC_FALSE; 2760 B->same_nonzero = PETSC_FALSE; 2761 2762 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 2763 2764 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 2765 "MatSeqAIJSetColumnIndices_SeqAIJ", 2766 MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 2767 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 2768 "MatStoreValues_SeqAIJ", 2769 MatStoreValues_SeqAIJ);CHKERRQ(ierr); 2770 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 2771 "MatRetrieveValues_SeqAIJ", 2772 MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 2773 #if !defined(PETSC_USE_64BIT_INT) 2774 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqsbaij_C", 2775 "MatConvert_SeqAIJ_SeqSBAIJ", 2776 MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 2777 #endif 2778 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqbaij_C", 2779 "MatConvert_SeqAIJ_SeqBAIJ", 2780 MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 2781 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsTranspose_C", 2782 "MatIsTranspose_SeqAIJ", 2783 MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 2784 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocation_C", 2785 "MatSeqAIJSetPreallocation_SeqAIJ", 2786 MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 2787 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatReorderForNonzeroDiagonal_C", 2788 "MatReorderForNonzeroDiagonal_SeqAIJ", 2789 MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 2790 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatAdjustForInodes_C", 2791 "MatAdjustForInodes_SeqAIJ", 2792 MatAdjustForInodes_SeqAIJ);CHKERRQ(ierr); 2793 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJGetInodeSizes_C", 2794 "MatSeqAIJGetInodeSizes_SeqAIJ", 2795 MatSeqAIJGetInodeSizes_SeqAIJ);CHKERRQ(ierr); 2796 PetscFunctionReturn(0); 2797 } 2798 EXTERN_C_END 2799 2800 #undef __FUNCT__ 2801 #define __FUNCT__ "MatDuplicate_SeqAIJ" 2802 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 2803 { 2804 Mat C; 2805 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 2806 PetscErrorCode ierr; 2807 PetscInt i,m = A->m; 2808 size_t len; 2809 2810 PetscFunctionBegin; 2811 *B = 0; 2812 ierr = MatCreate(A->comm,A->m,A->n,A->m,A->n,&C);CHKERRQ(ierr); 2813 ierr = MatSetType(C,A->type_name);CHKERRQ(ierr); 2814 ierr = PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 2815 2816 c = (Mat_SeqAIJ*)C->data; 2817 2818 C->factor = A->factor; 2819 C->lupivotthreshold = A->lupivotthreshold; 2820 2821 c->row = 0; 2822 c->col = 0; 2823 c->icol = 0; 2824 c->reallocs = 0; 2825 c->lu_shift = PETSC_FALSE; 2826 2827 C->assembled = PETSC_TRUE; 2828 2829 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&c->imax);CHKERRQ(ierr); 2830 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&c->ilen);CHKERRQ(ierr); 2831 for (i=0; i<m; i++) { 2832 c->imax[i] = a->imax[i]; 2833 c->ilen[i] = a->ilen[i]; 2834 } 2835 2836 /* allocate the matrix space */ 2837 c->singlemalloc = PETSC_TRUE; 2838 len = ((size_t) (m+1))*sizeof(PetscInt)+(a->i[m])*(sizeof(PetscScalar)+sizeof(PetscInt)); 2839 ierr = PetscMalloc(len,&c->a);CHKERRQ(ierr); 2840 c->j = (PetscInt*)(c->a + a->i[m] ); 2841 c->i = c->j + a->i[m]; 2842 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 2843 if (m > 0) { 2844 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(PetscInt));CHKERRQ(ierr); 2845 if (cpvalues == MAT_COPY_VALUES) { 2846 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 2847 } else { 2848 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 2849 } 2850 } 2851 2852 PetscLogObjectMemory(C,len+2*(m+1)*sizeof(PetscInt)+sizeof(struct _p_Mat)+sizeof(Mat_SeqAIJ)); 2853 c->sorted = a->sorted; 2854 c->ignorezeroentries = a->ignorezeroentries; 2855 c->roworiented = a->roworiented; 2856 c->nonew = a->nonew; 2857 if (a->diag) { 2858 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 2859 PetscLogObjectMemory(C,(m+1)*sizeof(PetscInt)); 2860 for (i=0; i<m; i++) { 2861 c->diag[i] = a->diag[i]; 2862 } 2863 } else c->diag = 0; 2864 c->solve_work = 0; 2865 c->inode.use = a->inode.use; 2866 c->inode.limit = a->inode.limit; 2867 c->inode.max_limit = a->inode.max_limit; 2868 if (a->inode.size){ 2869 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&c->inode.size);CHKERRQ(ierr); 2870 c->inode.node_count = a->inode.node_count; 2871 ierr = PetscMemcpy(c->inode.size,a->inode.size,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 2872 } else { 2873 c->inode.size = 0; 2874 c->inode.node_count = 0; 2875 } 2876 c->saved_values = 0; 2877 c->idiag = 0; 2878 c->ilu_preserve_row_sums = a->ilu_preserve_row_sums; 2879 c->ssor = 0; 2880 c->keepzeroedrows = a->keepzeroedrows; 2881 c->freedata = PETSC_TRUE; 2882 c->xtoy = 0; 2883 c->XtoY = 0; 2884 2885 c->nz = a->nz; 2886 c->maxnz = a->maxnz; 2887 C->preallocated = PETSC_TRUE; 2888 2889 c->compressedrow.use = a->compressedrow.use; 2890 c->compressedrow.nrows = a->compressedrow.nrows; 2891 c->compressedrow.checked = a->compressedrow.checked; 2892 if ( a->compressedrow.checked && a->compressedrow.use){ 2893 i = a->compressedrow.nrows; 2894 ierr = PetscMalloc((2*i+1)*sizeof(PetscInt),&c->compressedrow.i);CHKERRQ(ierr); 2895 c->compressedrow.rindex = c->compressedrow.i + i + 1; 2896 ierr = PetscMemcpy(c->compressedrow.i,a->compressedrow.i,(i+1)*sizeof(PetscInt));CHKERRQ(ierr); 2897 ierr = PetscMemcpy(c->compressedrow.rindex,a->compressedrow.rindex,i*sizeof(PetscInt));CHKERRQ(ierr); 2898 } else { 2899 c->compressedrow.use = PETSC_FALSE; 2900 c->compressedrow.i = PETSC_NULL; 2901 c->compressedrow.rindex = PETSC_NULL; 2902 } 2903 C->same_nonzero = A->same_nonzero; 2904 2905 *B = C; 2906 ierr = PetscFListDuplicate(A->qlist,&C->qlist);CHKERRQ(ierr); 2907 PetscFunctionReturn(0); 2908 } 2909 2910 #undef __FUNCT__ 2911 #define __FUNCT__ "MatLoad_SeqAIJ" 2912 PetscErrorCode MatLoad_SeqAIJ(PetscViewer viewer,const MatType type,Mat *A) 2913 { 2914 Mat_SeqAIJ *a; 2915 Mat B; 2916 PetscErrorCode ierr; 2917 PetscInt i,nz,header[4],*rowlengths = 0,M,N; 2918 int fd; 2919 PetscMPIInt size; 2920 MPI_Comm comm; 2921 2922 PetscFunctionBegin; 2923 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 2924 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2925 if (size > 1) SETERRQ(PETSC_ERR_ARG_SIZ,"view must have one processor"); 2926 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 2927 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 2928 if (header[0] != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 2929 M = header[1]; N = header[2]; nz = header[3]; 2930 2931 if (nz < 0) { 2932 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 2933 } 2934 2935 /* read in row lengths */ 2936 ierr = PetscMalloc(M*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 2937 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 2938 2939 /* create our matrix */ 2940 ierr = MatCreate(comm,PETSC_DECIDE,PETSC_DECIDE,M,N,&B);CHKERRQ(ierr); 2941 ierr = MatSetType(B,type);CHKERRQ(ierr); 2942 ierr = MatSeqAIJSetPreallocation(B,0,rowlengths);CHKERRQ(ierr); 2943 a = (Mat_SeqAIJ*)B->data; 2944 2945 /* read in column indices and adjust for Fortran indexing*/ 2946 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 2947 2948 /* read in nonzero values */ 2949 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 2950 2951 /* set matrix "i" values */ 2952 a->i[0] = 0; 2953 for (i=1; i<= M; i++) { 2954 a->i[i] = a->i[i-1] + rowlengths[i-1]; 2955 a->ilen[i-1] = rowlengths[i-1]; 2956 } 2957 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 2958 2959 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2960 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2961 *A = B; 2962 PetscFunctionReturn(0); 2963 } 2964 2965 #undef __FUNCT__ 2966 #define __FUNCT__ "MatEqual_SeqAIJ" 2967 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg) 2968 { 2969 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data,*b = (Mat_SeqAIJ *)B->data; 2970 PetscErrorCode ierr; 2971 2972 PetscFunctionBegin; 2973 /* If the matrix dimensions are not equal,or no of nonzeros */ 2974 if ((A->m != B->m) || (A->n != B->n) ||(a->nz != b->nz)) { 2975 *flg = PETSC_FALSE; 2976 PetscFunctionReturn(0); 2977 } 2978 2979 /* if the a->i are the same */ 2980 ierr = PetscMemcmp(a->i,b->i,(A->m+1)*sizeof(PetscInt),flg);CHKERRQ(ierr); 2981 if (!*flg) PetscFunctionReturn(0); 2982 2983 /* if a->j are the same */ 2984 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),flg);CHKERRQ(ierr); 2985 if (!*flg) PetscFunctionReturn(0); 2986 2987 /* if a->a are the same */ 2988 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 2989 2990 PetscFunctionReturn(0); 2991 2992 } 2993 2994 #undef __FUNCT__ 2995 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 2996 /*@C 2997 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 2998 provided by the user. 2999 3000 Coolective on MPI_Comm 3001 3002 Input Parameters: 3003 + comm - must be an MPI communicator of size 1 3004 . m - number of rows 3005 . n - number of columns 3006 . i - row indices 3007 . j - column indices 3008 - a - matrix values 3009 3010 Output Parameter: 3011 . mat - the matrix 3012 3013 Level: intermediate 3014 3015 Notes: 3016 The i, j, and a arrays are not copied by this routine, the user must free these arrays 3017 once the matrix is destroyed 3018 3019 You cannot set new nonzero locations into this matrix, that will generate an error. 3020 3021 The i and j indices are 0 based 3022 3023 .seealso: MatCreate(), MatCreateMPIAIJ(), MatCreateSeqAIJ() 3024 3025 @*/ 3026 PetscErrorCode MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt* i,PetscInt*j,PetscScalar *a,Mat *mat) 3027 { 3028 PetscErrorCode ierr; 3029 PetscInt ii; 3030 Mat_SeqAIJ *aij; 3031 3032 PetscFunctionBegin; 3033 ierr = MatCreate(comm,m,n,m,n,mat);CHKERRQ(ierr); 3034 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 3035 ierr = MatSeqAIJSetPreallocation(*mat,SKIP_ALLOCATION,0);CHKERRQ(ierr); 3036 aij = (Mat_SeqAIJ*)(*mat)->data; 3037 3038 if (i[0] != 0) { 3039 SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 3040 } 3041 aij->i = i; 3042 aij->j = j; 3043 aij->a = a; 3044 aij->singlemalloc = PETSC_FALSE; 3045 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 3046 aij->freedata = PETSC_FALSE; 3047 3048 for (ii=0; ii<m; ii++) { 3049 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 3050 #if defined(PETSC_USE_DEBUG) 3051 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]); 3052 #endif 3053 } 3054 #if defined(PETSC_USE_DEBUG) 3055 for (ii=0; ii<aij->i[m]; ii++) { 3056 if (j[ii] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 3057 if (j[ii] > n - 1) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column index to large at location = %d index = %d",ii,j[ii]); 3058 } 3059 #endif 3060 3061 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3062 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3063 PetscFunctionReturn(0); 3064 } 3065 3066 #undef __FUNCT__ 3067 #define __FUNCT__ "MatSetColoring_SeqAIJ" 3068 PetscErrorCode MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 3069 { 3070 PetscErrorCode ierr; 3071 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3072 3073 PetscFunctionBegin; 3074 if (coloring->ctype == IS_COLORING_LOCAL) { 3075 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 3076 a->coloring = coloring; 3077 } else if (coloring->ctype == IS_COLORING_GHOSTED) { 3078 PetscInt i,*larray; 3079 ISColoring ocoloring; 3080 ISColoringValue *colors; 3081 3082 /* set coloring for diagonal portion */ 3083 ierr = PetscMalloc((A->n+1)*sizeof(PetscInt),&larray);CHKERRQ(ierr); 3084 for (i=0; i<A->n; i++) { 3085 larray[i] = i; 3086 } 3087 ierr = ISGlobalToLocalMappingApply(A->mapping,IS_GTOLM_MASK,A->n,larray,PETSC_NULL,larray);CHKERRQ(ierr); 3088 ierr = PetscMalloc((A->n+1)*sizeof(ISColoringValue),&colors);CHKERRQ(ierr); 3089 for (i=0; i<A->n; i++) { 3090 colors[i] = coloring->colors[larray[i]]; 3091 } 3092 ierr = PetscFree(larray);CHKERRQ(ierr); 3093 ierr = ISColoringCreate(PETSC_COMM_SELF,A->n,colors,&ocoloring);CHKERRQ(ierr); 3094 a->coloring = ocoloring; 3095 } 3096 PetscFunctionReturn(0); 3097 } 3098 3099 #if defined(PETSC_HAVE_ADIC) 3100 EXTERN_C_BEGIN 3101 #include "adic/ad_utils.h" 3102 EXTERN_C_END 3103 3104 #undef __FUNCT__ 3105 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 3106 PetscErrorCode MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 3107 { 3108 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3109 PetscInt m = A->m,*ii = a->i,*jj = a->j,nz,i,j,nlen; 3110 PetscScalar *v = a->a,*values = ((PetscScalar*)advalues)+1; 3111 ISColoringValue *color; 3112 3113 PetscFunctionBegin; 3114 if (!a->coloring) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 3115 nlen = PetscADGetDerivTypeSize()/sizeof(PetscScalar); 3116 color = a->coloring->colors; 3117 /* loop over rows */ 3118 for (i=0; i<m; i++) { 3119 nz = ii[i+1] - ii[i]; 3120 /* loop over columns putting computed value into matrix */ 3121 for (j=0; j<nz; j++) { 3122 *v++ = values[color[*jj++]]; 3123 } 3124 values += nlen; /* jump to next row of derivatives */ 3125 } 3126 PetscFunctionReturn(0); 3127 } 3128 #endif 3129 3130 #undef __FUNCT__ 3131 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 3132 PetscErrorCode MatSetValuesAdifor_SeqAIJ(Mat A,PetscInt nl,void *advalues) 3133 { 3134 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3135 PetscInt m = A->m,*ii = a->i,*jj = a->j,nz,i,j; 3136 PetscScalar *v = a->a,*values = (PetscScalar *)advalues; 3137 ISColoringValue *color; 3138 3139 PetscFunctionBegin; 3140 if (!a->coloring) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 3141 color = a->coloring->colors; 3142 /* loop over rows */ 3143 for (i=0; i<m; i++) { 3144 nz = ii[i+1] - ii[i]; 3145 /* loop over columns putting computed value into matrix */ 3146 for (j=0; j<nz; j++) { 3147 *v++ = values[color[*jj++]]; 3148 } 3149 values += nl; /* jump to next row of derivatives */ 3150 } 3151 PetscFunctionReturn(0); 3152 } 3153 3154