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