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