1 #define PETSCMAT_DLL 2 3 4 /* 5 Defines the basic matrix operations for the AIJ (compressed row) 6 matrix storage format. 7 */ 8 9 10 #include "../src/mat/impls/aij/seq/aij.h" /*I "petscmat.h" I*/ 11 #include "petscbt.h" 12 13 #undef __FUNCT__ 14 #define __FUNCT__ "MatDiagonalSet_SeqAIJ" 15 PetscErrorCode PETSCMAT_DLLEXPORT MatDiagonalSet_SeqAIJ(Mat Y,Vec D,InsertMode is) 16 { 17 PetscErrorCode ierr; 18 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) Y->data; 19 PetscInt i,*diag, m = Y->rmap->n; 20 MatScalar *aa = aij->a; 21 PetscScalar *v; 22 PetscTruth missing; 23 24 PetscFunctionBegin; 25 if (Y->assembled) { 26 ierr = MatMissingDiagonal_SeqAIJ(Y,&missing,PETSC_NULL);CHKERRQ(ierr); 27 if (!missing) { 28 diag = aij->diag; 29 ierr = VecGetArray(D,&v);CHKERRQ(ierr); 30 if (is == INSERT_VALUES) { 31 for (i=0; i<m; i++) { 32 aa[diag[i]] = v[i]; 33 } 34 } else { 35 for (i=0; i<m; i++) { 36 aa[diag[i]] += v[i]; 37 } 38 } 39 ierr = VecRestoreArray(D,&v);CHKERRQ(ierr); 40 PetscFunctionReturn(0); 41 } 42 } 43 ierr = MatDiagonalSet_Default(Y,D,is);CHKERRQ(ierr); 44 PetscFunctionReturn(0); 45 } 46 47 #undef __FUNCT__ 48 #define __FUNCT__ "MatGetRowIJ_SeqAIJ" 49 PetscErrorCode MatGetRowIJ_SeqAIJ(Mat A,PetscInt oshift,PetscTruth symmetric,PetscTruth inodecompressed,PetscInt *m,PetscInt *ia[],PetscInt *ja[],PetscTruth *done) 50 { 51 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 52 PetscErrorCode ierr; 53 PetscInt i,ishift; 54 55 PetscFunctionBegin; 56 *m = A->rmap->n; 57 if (!ia) PetscFunctionReturn(0); 58 ishift = 0; 59 if (symmetric && !A->structurally_symmetric) { 60 ierr = MatToSymmetricIJ_SeqAIJ(A->rmap->n,a->i,a->j,ishift,oshift,ia,ja);CHKERRQ(ierr); 61 } else if (oshift == 1) { 62 PetscInt nz = a->i[A->rmap->n]; 63 /* malloc space and add 1 to i and j indices */ 64 ierr = PetscMalloc((A->rmap->n+1)*sizeof(PetscInt),ia);CHKERRQ(ierr); 65 for (i=0; i<A->rmap->n+1; i++) (*ia)[i] = a->i[i] + 1; 66 if (ja) { 67 ierr = PetscMalloc((nz+1)*sizeof(PetscInt),ja);CHKERRQ(ierr); 68 for (i=0; i<nz; i++) (*ja)[i] = a->j[i] + 1; 69 } 70 } else { 71 *ia = a->i; 72 if (ja) *ja = a->j; 73 } 74 PetscFunctionReturn(0); 75 } 76 77 #undef __FUNCT__ 78 #define __FUNCT__ "MatRestoreRowIJ_SeqAIJ" 79 PetscErrorCode MatRestoreRowIJ_SeqAIJ(Mat A,PetscInt oshift,PetscTruth symmetric,PetscTruth inodecompressed,PetscInt *n,PetscInt *ia[],PetscInt *ja[],PetscTruth *done) 80 { 81 PetscErrorCode ierr; 82 83 PetscFunctionBegin; 84 if (!ia) PetscFunctionReturn(0); 85 if ((symmetric && !A->structurally_symmetric) || oshift == 1) { 86 ierr = PetscFree(*ia);CHKERRQ(ierr); 87 if (ja) {ierr = PetscFree(*ja);CHKERRQ(ierr);} 88 } 89 PetscFunctionReturn(0); 90 } 91 92 #undef __FUNCT__ 93 #define __FUNCT__ "MatGetColumnIJ_SeqAIJ" 94 PetscErrorCode MatGetColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscTruth symmetric,PetscTruth inodecompressed,PetscInt *nn,PetscInt *ia[],PetscInt *ja[],PetscTruth *done) 95 { 96 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 97 PetscErrorCode ierr; 98 PetscInt i,*collengths,*cia,*cja,n = A->cmap->n,m = A->rmap->n; 99 PetscInt nz = a->i[m],row,*jj,mr,col; 100 101 PetscFunctionBegin; 102 *nn = n; 103 if (!ia) PetscFunctionReturn(0); 104 if (symmetric) { 105 ierr = MatToSymmetricIJ_SeqAIJ(A->rmap->n,a->i,a->j,0,oshift,ia,ja);CHKERRQ(ierr); 106 } else { 107 ierr = PetscMalloc((n+1)*sizeof(PetscInt),&collengths);CHKERRQ(ierr); 108 ierr = PetscMemzero(collengths,n*sizeof(PetscInt));CHKERRQ(ierr); 109 ierr = PetscMalloc((n+1)*sizeof(PetscInt),&cia);CHKERRQ(ierr); 110 ierr = PetscMalloc((nz+1)*sizeof(PetscInt),&cja);CHKERRQ(ierr); 111 jj = a->j; 112 for (i=0; i<nz; i++) { 113 collengths[jj[i]]++; 114 } 115 cia[0] = oshift; 116 for (i=0; i<n; i++) { 117 cia[i+1] = cia[i] + collengths[i]; 118 } 119 ierr = PetscMemzero(collengths,n*sizeof(PetscInt));CHKERRQ(ierr); 120 jj = a->j; 121 for (row=0; row<m; row++) { 122 mr = a->i[row+1] - a->i[row]; 123 for (i=0; i<mr; i++) { 124 col = *jj++; 125 cja[cia[col] + collengths[col]++ - oshift] = row + oshift; 126 } 127 } 128 ierr = PetscFree(collengths);CHKERRQ(ierr); 129 *ia = cia; *ja = cja; 130 } 131 PetscFunctionReturn(0); 132 } 133 134 #undef __FUNCT__ 135 #define __FUNCT__ "MatRestoreColumnIJ_SeqAIJ" 136 PetscErrorCode MatRestoreColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscTruth symmetric,PetscTruth inodecompressed,PetscInt *n,PetscInt *ia[],PetscInt *ja[],PetscTruth *done) 137 { 138 PetscErrorCode ierr; 139 140 PetscFunctionBegin; 141 if (!ia) PetscFunctionReturn(0); 142 143 ierr = PetscFree(*ia);CHKERRQ(ierr); 144 ierr = PetscFree(*ja);CHKERRQ(ierr); 145 146 PetscFunctionReturn(0); 147 } 148 149 #undef __FUNCT__ 150 #define __FUNCT__ "MatSetValuesRow_SeqAIJ" 151 PetscErrorCode MatSetValuesRow_SeqAIJ(Mat A,PetscInt row,const PetscScalar v[]) 152 { 153 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 154 PetscInt *ai = a->i; 155 PetscErrorCode ierr; 156 157 PetscFunctionBegin; 158 ierr = PetscMemcpy(a->a+ai[row],v,(ai[row+1]-ai[row])*sizeof(PetscScalar));CHKERRQ(ierr); 159 PetscFunctionReturn(0); 160 } 161 162 #define CHUNKSIZE 15 163 164 #undef __FUNCT__ 165 #define __FUNCT__ "MatSetValues_SeqAIJ" 166 PetscErrorCode MatSetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 167 { 168 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 169 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 170 PetscInt *imax = a->imax,*ai = a->i,*ailen = a->ilen; 171 PetscErrorCode ierr; 172 PetscInt *aj = a->j,nonew = a->nonew,lastcol = -1; 173 MatScalar *ap,value,*aa = a->a; 174 PetscTruth ignorezeroentries = a->ignorezeroentries; 175 PetscTruth roworiented = a->roworiented; 176 177 PetscFunctionBegin; 178 if (v) PetscValidScalarPointer(v,6); 179 for (k=0; k<m; k++) { /* loop over added rows */ 180 row = im[k]; 181 if (row < 0) continue; 182 #if defined(PETSC_USE_DEBUG) 183 if (row >= A->rmap->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->rmap->n-1); 184 #endif 185 rp = aj + ai[row]; ap = aa + ai[row]; 186 rmax = imax[row]; nrow = ailen[row]; 187 low = 0; 188 high = nrow; 189 for (l=0; l<n; l++) { /* loop over added columns */ 190 if (in[l] < 0) continue; 191 #if defined(PETSC_USE_DEBUG) 192 if (in[l] >= A->cmap->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->cmap->n-1); 193 #endif 194 col = in[l]; 195 if (v) { 196 if (roworiented) { 197 value = v[l + k*n]; 198 } else { 199 value = v[k + l*m]; 200 } 201 } else { 202 value = 0.; 203 } 204 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 205 206 if (col <= lastcol) low = 0; else high = nrow; 207 lastcol = col; 208 while (high-low > 5) { 209 t = (low+high)/2; 210 if (rp[t] > col) high = t; 211 else low = t; 212 } 213 for (i=low; i<high; i++) { 214 if (rp[i] > col) break; 215 if (rp[i] == col) { 216 if (is == ADD_VALUES) ap[i] += value; 217 else ap[i] = value; 218 low = i + 1; 219 goto noinsert; 220 } 221 } 222 if (value == 0.0 && ignorezeroentries) goto noinsert; 223 if (nonew == 1) goto noinsert; 224 if (nonew == -1) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at (%D,%D) in the matrix",row,col); 225 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 226 N = nrow++ - 1; a->nz++; high++; 227 /* shift up all the later entries in this row */ 228 for (ii=N; ii>=i; ii--) { 229 rp[ii+1] = rp[ii]; 230 ap[ii+1] = ap[ii]; 231 } 232 rp[i] = col; 233 ap[i] = value; 234 low = i + 1; 235 noinsert:; 236 } 237 ailen[row] = nrow; 238 } 239 A->same_nonzero = PETSC_FALSE; 240 PetscFunctionReturn(0); 241 } 242 243 244 #undef __FUNCT__ 245 #define __FUNCT__ "MatGetValues_SeqAIJ" 246 PetscErrorCode MatGetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],PetscScalar v[]) 247 { 248 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 249 PetscInt *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j; 250 PetscInt *ai = a->i,*ailen = a->ilen; 251 MatScalar *ap,*aa = a->a; 252 253 PetscFunctionBegin; 254 for (k=0; k<m; k++) { /* loop over rows */ 255 row = im[k]; 256 if (row < 0) {v += n; continue;} /* SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Negative row: %D",row); */ 257 if (row >= A->rmap->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->rmap->n-1); 258 rp = aj + ai[row]; ap = aa + ai[row]; 259 nrow = ailen[row]; 260 for (l=0; l<n; l++) { /* loop over columns */ 261 if (in[l] < 0) {v++; continue;} /* SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Negative column: %D",in[l]); */ 262 if (in[l] >= A->cmap->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->cmap->n-1); 263 col = in[l] ; 264 high = nrow; low = 0; /* assume unsorted */ 265 while (high-low > 5) { 266 t = (low+high)/2; 267 if (rp[t] > col) high = t; 268 else low = t; 269 } 270 for (i=low; i<high; i++) { 271 if (rp[i] > col) break; 272 if (rp[i] == col) { 273 *v++ = ap[i]; 274 goto finished; 275 } 276 } 277 *v++ = 0.0; 278 finished:; 279 } 280 } 281 PetscFunctionReturn(0); 282 } 283 284 285 #undef __FUNCT__ 286 #define __FUNCT__ "MatView_SeqAIJ_Binary" 287 PetscErrorCode MatView_SeqAIJ_Binary(Mat A,PetscViewer viewer) 288 { 289 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 290 PetscErrorCode ierr; 291 PetscInt i,*col_lens; 292 int fd; 293 294 PetscFunctionBegin; 295 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 296 ierr = PetscMalloc((4+A->rmap->n)*sizeof(PetscInt),&col_lens);CHKERRQ(ierr); 297 col_lens[0] = MAT_FILE_COOKIE; 298 col_lens[1] = A->rmap->n; 299 col_lens[2] = A->cmap->n; 300 col_lens[3] = a->nz; 301 302 /* store lengths of each row and write (including header) to file */ 303 for (i=0; i<A->rmap->n; i++) { 304 col_lens[4+i] = a->i[i+1] - a->i[i]; 305 } 306 ierr = PetscBinaryWrite(fd,col_lens,4+A->rmap->n,PETSC_INT,PETSC_TRUE);CHKERRQ(ierr); 307 ierr = PetscFree(col_lens);CHKERRQ(ierr); 308 309 /* store column indices (zero start index) */ 310 ierr = PetscBinaryWrite(fd,a->j,a->nz,PETSC_INT,PETSC_FALSE);CHKERRQ(ierr); 311 312 /* store nonzero values */ 313 ierr = PetscBinaryWrite(fd,a->a,a->nz,PETSC_SCALAR,PETSC_FALSE);CHKERRQ(ierr); 314 PetscFunctionReturn(0); 315 } 316 317 EXTERN PetscErrorCode MatSeqAIJFactorInfo_Matlab(Mat,PetscViewer); 318 319 #undef __FUNCT__ 320 #define __FUNCT__ "MatView_SeqAIJ_ASCII" 321 PetscErrorCode MatView_SeqAIJ_ASCII(Mat A,PetscViewer viewer) 322 { 323 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 324 PetscErrorCode ierr; 325 PetscInt i,j,m = A->rmap->n,shift=0; 326 const char *name; 327 PetscViewerFormat format; 328 329 PetscFunctionBegin; 330 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 331 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 332 if (format == PETSC_VIEWER_ASCII_MATLAB) { 333 PetscInt nofinalvalue = 0; 334 if ((a->i[m] == a->i[m-1]) || (a->j[a->nz-1] != A->cmap->n-!shift)) { 335 nofinalvalue = 1; 336 } 337 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 338 ierr = PetscViewerASCIIPrintf(viewer,"%% Size = %D %D \n",m,A->cmap->n);CHKERRQ(ierr); 339 ierr = PetscViewerASCIIPrintf(viewer,"%% Nonzeros = %D \n",a->nz);CHKERRQ(ierr); 340 ierr = PetscViewerASCIIPrintf(viewer,"zzz = zeros(%D,3);\n",a->nz+nofinalvalue);CHKERRQ(ierr); 341 ierr = PetscViewerASCIIPrintf(viewer,"zzz = [\n");CHKERRQ(ierr); 342 343 for (i=0; i<m; i++) { 344 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 345 #if defined(PETSC_USE_COMPLEX) 346 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); 347 #else 348 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %18.16e\n",i+1,a->j[j]+!shift,a->a[j]);CHKERRQ(ierr); 349 #endif 350 } 351 } 352 if (nofinalvalue) { 353 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %18.16e\n",m,A->cmap->n,0.0);CHKERRQ(ierr); 354 } 355 ierr = PetscViewerASCIIPrintf(viewer,"];\n %s = spconvert(zzz);\n",name);CHKERRQ(ierr); 356 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 357 } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO) { 358 PetscFunctionReturn(0); 359 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 360 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 361 for (i=0; i<m; i++) { 362 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 363 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 364 #if defined(PETSC_USE_COMPLEX) 365 if (PetscImaginaryPart(a->a[j]) > 0.0 && PetscRealPart(a->a[j]) != 0.0) { 366 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 367 } else if (PetscImaginaryPart(a->a[j]) < 0.0 && PetscRealPart(a->a[j]) != 0.0) { 368 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 369 } else if (PetscRealPart(a->a[j]) != 0.0) { 370 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 371 } 372 #else 373 if (a->a[j] != 0.0) {ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr);} 374 #endif 375 } 376 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 377 } 378 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 379 } else if (format == PETSC_VIEWER_ASCII_SYMMODU) { 380 PetscInt nzd=0,fshift=1,*sptr; 381 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 382 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&sptr);CHKERRQ(ierr); 383 for (i=0; i<m; i++) { 384 sptr[i] = nzd+1; 385 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 386 if (a->j[j] >= i) { 387 #if defined(PETSC_USE_COMPLEX) 388 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) nzd++; 389 #else 390 if (a->a[j] != 0.0) nzd++; 391 #endif 392 } 393 } 394 } 395 sptr[m] = nzd+1; 396 ierr = PetscViewerASCIIPrintf(viewer," %D %D\n\n",m,nzd);CHKERRQ(ierr); 397 for (i=0; i<m+1; i+=6) { 398 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);} 399 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);} 400 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);} 401 else if (i+1<m) {ierr = PetscViewerASCIIPrintf(viewer," %D %D %D\n",sptr[i],sptr[i+1],sptr[i+2]);CHKERRQ(ierr);} 402 else if (i<m) {ierr = PetscViewerASCIIPrintf(viewer," %D %D\n",sptr[i],sptr[i+1]);CHKERRQ(ierr);} 403 else {ierr = PetscViewerASCIIPrintf(viewer," %D\n",sptr[i]);CHKERRQ(ierr);} 404 } 405 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 406 ierr = PetscFree(sptr);CHKERRQ(ierr); 407 for (i=0; i<m; i++) { 408 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 409 if (a->j[j] >= i) {ierr = PetscViewerASCIIPrintf(viewer," %D ",a->j[j]+fshift);CHKERRQ(ierr);} 410 } 411 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 412 } 413 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 414 for (i=0; i<m; i++) { 415 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 416 if (a->j[j] >= i) { 417 #if defined(PETSC_USE_COMPLEX) 418 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) { 419 ierr = PetscViewerASCIIPrintf(viewer," %18.16e %18.16e ",PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 420 } 421 #else 422 if (a->a[j] != 0.0) {ierr = PetscViewerASCIIPrintf(viewer," %18.16e ",a->a[j]);CHKERRQ(ierr);} 423 #endif 424 } 425 } 426 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 427 } 428 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 429 } else if (format == PETSC_VIEWER_ASCII_DENSE) { 430 PetscInt cnt = 0,jcnt; 431 PetscScalar value; 432 433 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 434 for (i=0; i<m; i++) { 435 jcnt = 0; 436 for (j=0; j<A->cmap->n; j++) { 437 if (jcnt < a->i[i+1]-a->i[i] && j == a->j[cnt]) { 438 value = a->a[cnt++]; 439 jcnt++; 440 } else { 441 value = 0.0; 442 } 443 #if defined(PETSC_USE_COMPLEX) 444 ierr = PetscViewerASCIIPrintf(viewer," %7.5e+%7.5e i ",PetscRealPart(value),PetscImaginaryPart(value));CHKERRQ(ierr); 445 #else 446 ierr = PetscViewerASCIIPrintf(viewer," %7.5e ",value);CHKERRQ(ierr); 447 #endif 448 } 449 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 450 } 451 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 452 } else if (format == PETSC_VIEWER_ASCII_MATRIXMARKET) { 453 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 454 #if defined(PETSC_USE_COMPLEX) 455 ierr = PetscViewerASCIIPrintf(viewer,"%%matrix complex general\n");CHKERRQ(ierr); 456 #else 457 ierr = PetscViewerASCIIPrintf(viewer,"%%matrix real general\n");CHKERRQ(ierr); 458 #endif 459 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %D\n", m, A->cmap->n, a->nz);CHKERRQ(ierr); 460 for (i=0; i<m; i++) { 461 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 462 #if defined(PETSC_USE_COMPLEX) 463 if (PetscImaginaryPart(a->a[j]) > 0.0) { 464 ierr = PetscViewerASCIIPrintf(viewer,"%D %D, %G %G\n", i+shift,a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 465 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 466 ierr = PetscViewerASCIIPrintf(viewer,"%D %D, %G -%G\n", i+shift,a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 467 } else { 468 ierr = PetscViewerASCIIPrintf(viewer,"%D %D, %G\n", i+shift,a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 469 } 470 #else 471 ierr = PetscViewerASCIIPrintf(viewer,"%D %D %G\n", i+shift, a->j[j]+shift, a->a[j]);CHKERRQ(ierr); 472 #endif 473 } 474 } 475 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 476 } else { 477 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 478 if (A->factor){ 479 for (i=0; i<m; i++) { 480 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 481 /* L part */ 482 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 483 #if defined(PETSC_USE_COMPLEX) 484 if (PetscImaginaryPart(a->a[j]) > 0.0) { 485 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 486 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 487 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 488 } else { 489 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 490 } 491 #else 492 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr); 493 #endif 494 } 495 /* diagonal */ 496 j = a->diag[i]; 497 #if defined(PETSC_USE_COMPLEX) 498 if (PetscImaginaryPart(a->a[j]) > 0.0) { 499 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 500 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 501 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 502 } else { 503 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 504 } 505 #else 506 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr); 507 #endif 508 509 /* U part */ 510 for (j=a->diag[i+1]+1+shift; j<a->diag[i]+shift; j++) { 511 #if defined(PETSC_USE_COMPLEX) 512 if (PetscImaginaryPart(a->a[j]) > 0.0) { 513 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 514 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 515 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 516 } else { 517 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 518 } 519 #else 520 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr); 521 #endif 522 } 523 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 524 } 525 } else { 526 for (i=0; i<m; i++) { 527 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 528 for (j=a->i[i]+shift; j<a->i[i+1]+shift; j++) { 529 #if defined(PETSC_USE_COMPLEX) 530 if (PetscImaginaryPart(a->a[j]) > 0.0) { 531 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 532 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 533 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i)",a->j[j]+shift,PetscRealPart(a->a[j]),-PetscImaginaryPart(a->a[j]));CHKERRQ(ierr); 534 } else { 535 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,PetscRealPart(a->a[j]));CHKERRQ(ierr); 536 } 537 #else 538 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[j]+shift,a->a[j]);CHKERRQ(ierr); 539 #endif 540 } 541 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 542 } 543 } 544 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 545 } 546 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 547 PetscFunctionReturn(0); 548 } 549 550 #undef __FUNCT__ 551 #define __FUNCT__ "MatView_SeqAIJ_Draw_Zoom" 552 PetscErrorCode MatView_SeqAIJ_Draw_Zoom(PetscDraw draw,void *Aa) 553 { 554 Mat A = (Mat) Aa; 555 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 556 PetscErrorCode ierr; 557 PetscInt i,j,m = A->rmap->n,color; 558 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r,maxv = 0.0; 559 PetscViewer viewer; 560 PetscViewerFormat format; 561 562 PetscFunctionBegin; 563 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 564 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 565 566 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 567 /* loop over matrix elements drawing boxes */ 568 569 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 570 /* Blue for negative, Cyan for zero and Red for positive */ 571 color = PETSC_DRAW_BLUE; 572 for (i=0; i<m; i++) { 573 y_l = m - i - 1.0; y_r = y_l + 1.0; 574 for (j=a->i[i]; j<a->i[i+1]; j++) { 575 x_l = a->j[j] ; x_r = x_l + 1.0; 576 #if defined(PETSC_USE_COMPLEX) 577 if (PetscRealPart(a->a[j]) >= 0.) continue; 578 #else 579 if (a->a[j] >= 0.) continue; 580 #endif 581 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 582 } 583 } 584 color = PETSC_DRAW_CYAN; 585 for (i=0; i<m; i++) { 586 y_l = m - i - 1.0; y_r = y_l + 1.0; 587 for (j=a->i[i]; j<a->i[i+1]; j++) { 588 x_l = a->j[j]; x_r = x_l + 1.0; 589 if (a->a[j] != 0.) continue; 590 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 591 } 592 } 593 color = PETSC_DRAW_RED; 594 for (i=0; i<m; i++) { 595 y_l = m - i - 1.0; y_r = y_l + 1.0; 596 for (j=a->i[i]; j<a->i[i+1]; j++) { 597 x_l = a->j[j]; x_r = x_l + 1.0; 598 #if defined(PETSC_USE_COMPLEX) 599 if (PetscRealPart(a->a[j]) <= 0.) continue; 600 #else 601 if (a->a[j] <= 0.) continue; 602 #endif 603 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 604 } 605 } 606 } else { 607 /* use contour shading to indicate magnitude of values */ 608 /* first determine max of all nonzero values */ 609 PetscInt nz = a->nz,count; 610 PetscDraw popup; 611 PetscReal scale; 612 613 for (i=0; i<nz; i++) { 614 if (PetscAbsScalar(a->a[i]) > maxv) maxv = PetscAbsScalar(a->a[i]); 615 } 616 scale = (245.0 - PETSC_DRAW_BASIC_COLORS)/maxv; 617 ierr = PetscDrawGetPopup(draw,&popup);CHKERRQ(ierr); 618 if (popup) {ierr = PetscDrawScalePopup(popup,0.0,maxv);CHKERRQ(ierr);} 619 count = 0; 620 for (i=0; i<m; i++) { 621 y_l = m - i - 1.0; y_r = y_l + 1.0; 622 for (j=a->i[i]; j<a->i[i+1]; j++) { 623 x_l = a->j[j]; x_r = x_l + 1.0; 624 color = PETSC_DRAW_BASIC_COLORS + (PetscInt)(scale*PetscAbsScalar(a->a[count])); 625 ierr = PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color);CHKERRQ(ierr); 626 count++; 627 } 628 } 629 } 630 PetscFunctionReturn(0); 631 } 632 633 #undef __FUNCT__ 634 #define __FUNCT__ "MatView_SeqAIJ_Draw" 635 PetscErrorCode MatView_SeqAIJ_Draw(Mat A,PetscViewer viewer) 636 { 637 PetscErrorCode ierr; 638 PetscDraw draw; 639 PetscReal xr,yr,xl,yl,h,w; 640 PetscTruth isnull; 641 642 PetscFunctionBegin; 643 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 644 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); 645 if (isnull) PetscFunctionReturn(0); 646 647 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 648 xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0; 649 xr += w; yr += h; xl = -w; yl = -h; 650 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 651 ierr = PetscDrawZoom(draw,MatView_SeqAIJ_Draw_Zoom,A);CHKERRQ(ierr); 652 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);CHKERRQ(ierr); 653 PetscFunctionReturn(0); 654 } 655 656 #undef __FUNCT__ 657 #define __FUNCT__ "MatView_SeqAIJ" 658 PetscErrorCode MatView_SeqAIJ(Mat A,PetscViewer viewer) 659 { 660 PetscErrorCode ierr; 661 PetscTruth iascii,isbinary,isdraw; 662 663 PetscFunctionBegin; 664 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 665 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_BINARY,&isbinary);CHKERRQ(ierr); 666 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_DRAW,&isdraw);CHKERRQ(ierr); 667 if (iascii) { 668 ierr = MatView_SeqAIJ_ASCII(A,viewer);CHKERRQ(ierr); 669 } else if (isbinary) { 670 ierr = MatView_SeqAIJ_Binary(A,viewer);CHKERRQ(ierr); 671 } else if (isdraw) { 672 ierr = MatView_SeqAIJ_Draw(A,viewer);CHKERRQ(ierr); 673 } else { 674 SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported by SeqAIJ matrices",((PetscObject)viewer)->type_name); 675 } 676 ierr = MatView_SeqAIJ_Inode(A,viewer);CHKERRQ(ierr); 677 PetscFunctionReturn(0); 678 } 679 680 #undef __FUNCT__ 681 #define __FUNCT__ "MatAssemblyEnd_SeqAIJ" 682 PetscErrorCode MatAssemblyEnd_SeqAIJ(Mat A,MatAssemblyType mode) 683 { 684 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 685 PetscErrorCode ierr; 686 PetscInt fshift = 0,i,j,*ai = a->i,*aj = a->j,*imax = a->imax; 687 PetscInt m = A->rmap->n,*ip,N,*ailen = a->ilen,rmax = 0; 688 MatScalar *aa = a->a,*ap; 689 PetscReal ratio=0.6; 690 691 PetscFunctionBegin; 692 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 693 694 if (m) rmax = ailen[0]; /* determine row with most nonzeros */ 695 for (i=1; i<m; i++) { 696 /* move each row back by the amount of empty slots (fshift) before it*/ 697 fshift += imax[i-1] - ailen[i-1]; 698 rmax = PetscMax(rmax,ailen[i]); 699 if (fshift) { 700 ip = aj + ai[i] ; 701 ap = aa + ai[i] ; 702 N = ailen[i]; 703 for (j=0; j<N; j++) { 704 ip[j-fshift] = ip[j]; 705 ap[j-fshift] = ap[j]; 706 } 707 } 708 ai[i] = ai[i-1] + ailen[i-1]; 709 } 710 if (m) { 711 fshift += imax[m-1] - ailen[m-1]; 712 ai[m] = ai[m-1] + ailen[m-1]; 713 } 714 /* reset ilen and imax for each row */ 715 for (i=0; i<m; i++) { 716 ailen[i] = imax[i] = ai[i+1] - ai[i]; 717 } 718 a->nz = ai[m]; 719 if (fshift && a->nounused == -1) { 720 SETERRQ3(PETSC_ERR_PLIB, "Unused space detected in matrix: %D X %D, %D unneeded", m, A->cmap->n, fshift); 721 } 722 723 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 724 ierr = PetscInfo4(A,"Matrix size: %D X %D; storage space: %D unneeded,%D used\n",m,A->cmap->n,fshift,a->nz);CHKERRQ(ierr); 725 ierr = PetscInfo1(A,"Number of mallocs during MatSetValues() is %D\n",a->reallocs);CHKERRQ(ierr); 726 ierr = PetscInfo1(A,"Maximum nonzeros in any row is %D\n",rmax);CHKERRQ(ierr); 727 728 a->reallocs = 0; 729 A->info.nz_unneeded = (double)fshift; 730 a->rmax = rmax; 731 732 /* check for zero rows. If found a large number of zero rows, use CompressedRow functions */ 733 ierr = Mat_CheckCompressedRow(A,&a->compressedrow,a->i,m,ratio);CHKERRQ(ierr); 734 A->same_nonzero = PETSC_TRUE; 735 736 ierr = MatAssemblyEnd_SeqAIJ_Inode(A,mode);CHKERRQ(ierr); 737 738 a->idiagvalid = PETSC_FALSE; 739 PetscFunctionReturn(0); 740 } 741 742 #undef __FUNCT__ 743 #define __FUNCT__ "MatRealPart_SeqAIJ" 744 PetscErrorCode MatRealPart_SeqAIJ(Mat A) 745 { 746 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 747 PetscInt i,nz = a->nz; 748 MatScalar *aa = a->a; 749 750 PetscFunctionBegin; 751 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 752 PetscFunctionReturn(0); 753 } 754 755 #undef __FUNCT__ 756 #define __FUNCT__ "MatImaginaryPart_SeqAIJ" 757 PetscErrorCode MatImaginaryPart_SeqAIJ(Mat A) 758 { 759 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 760 PetscInt i,nz = a->nz; 761 MatScalar *aa = a->a; 762 763 PetscFunctionBegin; 764 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 765 PetscFunctionReturn(0); 766 } 767 768 #undef __FUNCT__ 769 #define __FUNCT__ "MatZeroEntries_SeqAIJ" 770 PetscErrorCode MatZeroEntries_SeqAIJ(Mat A) 771 { 772 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 773 PetscErrorCode ierr; 774 775 PetscFunctionBegin; 776 ierr = PetscMemzero(a->a,(a->i[A->rmap->n])*sizeof(PetscScalar));CHKERRQ(ierr); 777 PetscFunctionReturn(0); 778 } 779 780 #undef __FUNCT__ 781 #define __FUNCT__ "MatDestroy_SeqAIJ" 782 PetscErrorCode MatDestroy_SeqAIJ(Mat A) 783 { 784 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 785 PetscErrorCode ierr; 786 787 PetscFunctionBegin; 788 #if defined(PETSC_USE_LOG) 789 PetscLogObjectState((PetscObject)A,"Rows=%D, Cols=%D, NZ=%D",A->rmap->n,A->cmap->n,a->nz); 790 #endif 791 ierr = MatSeqXAIJFreeAIJ(A,&a->a,&a->j,&a->i);CHKERRQ(ierr); 792 if (a->row) { 793 ierr = ISDestroy(a->row);CHKERRQ(ierr); 794 } 795 if (a->col) { 796 ierr = ISDestroy(a->col);CHKERRQ(ierr); 797 } 798 ierr = PetscFree(a->diag);CHKERRQ(ierr); 799 ierr = PetscFree2(a->imax,a->ilen);CHKERRQ(ierr); 800 ierr = PetscFree3(a->idiag,a->mdiag,a->ssor_work);CHKERRQ(ierr); 801 ierr = PetscFree(a->solve_work);CHKERRQ(ierr); 802 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} 803 ierr = PetscFree(a->saved_values);CHKERRQ(ierr); 804 if (a->coloring) {ierr = ISColoringDestroy(a->coloring);CHKERRQ(ierr);} 805 ierr = PetscFree(a->xtoy);CHKERRQ(ierr); 806 if (a->XtoY) {ierr = MatDestroy(a->XtoY);CHKERRQ(ierr);} 807 if (a->compressedrow.checked && a->compressedrow.use){ierr = PetscFree2(a->compressedrow.i,a->compressedrow.rindex);} 808 809 ierr = MatDestroy_SeqAIJ_Inode(A);CHKERRQ(ierr); 810 811 ierr = PetscFree(a);CHKERRQ(ierr); 812 813 ierr = PetscObjectChangeTypeName((PetscObject)A,0);CHKERRQ(ierr); 814 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatSeqAIJSetColumnIndices_C","",PETSC_NULL);CHKERRQ(ierr); 815 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatStoreValues_C","",PETSC_NULL);CHKERRQ(ierr); 816 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatRetrieveValues_C","",PETSC_NULL);CHKERRQ(ierr); 817 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatConvert_seqaij_seqsbaij_C","",PETSC_NULL);CHKERRQ(ierr); 818 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatConvert_seqaij_seqbaij_C","",PETSC_NULL);CHKERRQ(ierr); 819 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatConvert_seqaij_seqcsrperm_C","",PETSC_NULL);CHKERRQ(ierr); 820 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatIsTranspose_C","",PETSC_NULL);CHKERRQ(ierr); 821 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatSeqAIJSetPreallocation_C","",PETSC_NULL);CHKERRQ(ierr); 822 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatSeqAIJSetPreallocationCSR_C","",PETSC_NULL);CHKERRQ(ierr); 823 ierr = PetscObjectComposeFunctionDynamic((PetscObject)A,"MatReorderForNonzeroDiagonal_C","",PETSC_NULL);CHKERRQ(ierr); 824 PetscFunctionReturn(0); 825 } 826 827 #undef __FUNCT__ 828 #define __FUNCT__ "MatSetOption_SeqAIJ" 829 PetscErrorCode MatSetOption_SeqAIJ(Mat A,MatOption op,PetscTruth flg) 830 { 831 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 832 PetscErrorCode ierr; 833 834 PetscFunctionBegin; 835 switch (op) { 836 case MAT_ROW_ORIENTED: 837 a->roworiented = flg; 838 break; 839 case MAT_KEEP_NONZERO_PATTERN: 840 a->keepnonzeropattern = flg; 841 break; 842 case MAT_NEW_NONZERO_LOCATIONS: 843 a->nonew = (flg ? 0 : 1); 844 break; 845 case MAT_NEW_NONZERO_LOCATION_ERR: 846 a->nonew = (flg ? -1 : 0); 847 break; 848 case MAT_NEW_NONZERO_ALLOCATION_ERR: 849 a->nonew = (flg ? -2 : 0); 850 break; 851 case MAT_UNUSED_NONZERO_LOCATION_ERR: 852 a->nounused = (flg ? -1 : 0); 853 break; 854 case MAT_IGNORE_ZERO_ENTRIES: 855 a->ignorezeroentries = flg; 856 break; 857 case MAT_USE_COMPRESSEDROW: 858 a->compressedrow.use = flg; 859 break; 860 case MAT_SYMMETRIC: 861 case MAT_STRUCTURALLY_SYMMETRIC: 862 case MAT_HERMITIAN: 863 case MAT_SYMMETRY_ETERNAL: 864 case MAT_NEW_DIAGONALS: 865 case MAT_IGNORE_OFF_PROC_ENTRIES: 866 case MAT_USE_HASH_TABLE: 867 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 868 break; 869 case MAT_USE_INODES: 870 /* Not an error because MatSetOption_SeqAIJ_Inode handles this one */ 871 break; 872 default: 873 SETERRQ1(PETSC_ERR_SUP,"unknown option %d",op); 874 } 875 ierr = MatSetOption_SeqAIJ_Inode(A,op,flg);CHKERRQ(ierr); 876 PetscFunctionReturn(0); 877 } 878 879 #undef __FUNCT__ 880 #define __FUNCT__ "MatGetDiagonal_SeqAIJ" 881 PetscErrorCode MatGetDiagonal_SeqAIJ(Mat A,Vec v) 882 { 883 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 884 PetscErrorCode ierr; 885 PetscInt i,j,n; 886 PetscScalar *x,zero = 0.0; 887 888 PetscFunctionBegin; 889 ierr = VecSet(v,zero);CHKERRQ(ierr); 890 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 891 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 892 if (n != A->rmap->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 893 for (i=0; i<A->rmap->n; i++) { 894 for (j=a->i[i]; j<a->i[i+1]; j++) { 895 if (a->j[j] == i) { 896 x[i] = a->a[j]; 897 break; 898 } 899 } 900 } 901 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 902 PetscFunctionReturn(0); 903 } 904 905 #include "../src/mat/impls/aij/seq/ftn-kernels/fmult.h" 906 #undef __FUNCT__ 907 #define __FUNCT__ "MatMultTransposeAdd_SeqAIJ" 908 PetscErrorCode MatMultTransposeAdd_SeqAIJ(Mat A,Vec xx,Vec zz,Vec yy) 909 { 910 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 911 PetscScalar *x,*y; 912 PetscErrorCode ierr; 913 PetscInt m = A->rmap->n; 914 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 915 MatScalar *v; 916 PetscScalar alpha; 917 PetscInt n,i,j,*idx,*ii,*ridx=PETSC_NULL; 918 Mat_CompressedRow cprow = a->compressedrow; 919 PetscTruth usecprow = cprow.use; 920 #endif 921 922 PetscFunctionBegin; 923 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 924 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 925 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 926 927 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 928 fortranmulttransposeaddaij_(&m,x,a->i,a->j,a->a,y); 929 #else 930 if (usecprow){ 931 m = cprow.nrows; 932 ii = cprow.i; 933 ridx = cprow.rindex; 934 } else { 935 ii = a->i; 936 } 937 for (i=0; i<m; i++) { 938 idx = a->j + ii[i] ; 939 v = a->a + ii[i] ; 940 n = ii[i+1] - ii[i]; 941 if (usecprow){ 942 alpha = x[ridx[i]]; 943 } else { 944 alpha = x[i]; 945 } 946 for (j=0; j<n; j++) y[idx[j]] += alpha*v[j]; 947 } 948 #endif 949 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 950 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 951 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 952 PetscFunctionReturn(0); 953 } 954 955 #undef __FUNCT__ 956 #define __FUNCT__ "MatMultTranspose_SeqAIJ" 957 PetscErrorCode MatMultTranspose_SeqAIJ(Mat A,Vec xx,Vec yy) 958 { 959 PetscErrorCode ierr; 960 961 PetscFunctionBegin; 962 ierr = VecSet(yy,0.0);CHKERRQ(ierr); 963 ierr = MatMultTransposeAdd_SeqAIJ(A,xx,yy,yy);CHKERRQ(ierr); 964 PetscFunctionReturn(0); 965 } 966 967 #include "../src/mat/impls/aij/seq/ftn-kernels/fmult.h" 968 #undef __FUNCT__ 969 #define __FUNCT__ "MatMult_SeqAIJ" 970 PetscErrorCode MatMult_SeqAIJ(Mat A,Vec xx,Vec yy) 971 { 972 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 973 PetscScalar *y; 974 const PetscScalar *x; 975 const MatScalar *aa; 976 PetscErrorCode ierr; 977 PetscInt m=A->rmap->n; 978 const PetscInt *aj,*ii,*ridx=PETSC_NULL; 979 PetscInt n,i,nonzerorow=0; 980 PetscScalar sum; 981 PetscTruth usecprow=a->compressedrow.use; 982 983 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 984 #pragma disjoint(*x,*y,*aa) 985 #endif 986 987 PetscFunctionBegin; 988 ierr = VecGetArray(xx,(PetscScalar**)&x);CHKERRQ(ierr); 989 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 990 aj = a->j; 991 aa = a->a; 992 ii = a->i; 993 if (usecprow){ /* use compressed row format */ 994 m = a->compressedrow.nrows; 995 ii = a->compressedrow.i; 996 ridx = a->compressedrow.rindex; 997 for (i=0; i<m; i++){ 998 n = ii[i+1] - ii[i]; 999 aj = a->j + ii[i]; 1000 aa = a->a + ii[i]; 1001 sum = 0.0; 1002 nonzerorow += (n>0); 1003 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1004 /* for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; */ 1005 y[*ridx++] = sum; 1006 } 1007 } else { /* do not use compressed row format */ 1008 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 1009 fortranmultaij_(&m,x,ii,aj,aa,y); 1010 #else 1011 for (i=0; i<m; i++) { 1012 n = ii[i+1] - ii[i]; 1013 aj = a->j + ii[i]; 1014 aa = a->a + ii[i]; 1015 sum = 0.0; 1016 nonzerorow += (n>0); 1017 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1018 y[i] = sum; 1019 } 1020 #endif 1021 } 1022 ierr = PetscLogFlops(2.0*a->nz - nonzerorow);CHKERRQ(ierr); 1023 ierr = VecRestoreArray(xx,(PetscScalar**)&x);CHKERRQ(ierr); 1024 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1025 PetscFunctionReturn(0); 1026 } 1027 1028 #include "../src/mat/impls/aij/seq/ftn-kernels/fmultadd.h" 1029 #undef __FUNCT__ 1030 #define __FUNCT__ "MatMultAdd_SeqAIJ" 1031 PetscErrorCode MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 1032 { 1033 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1034 PetscScalar *x,*y,*z; 1035 const MatScalar *aa; 1036 PetscErrorCode ierr; 1037 PetscInt m = A->rmap->n,*aj,*ii; 1038 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 1039 PetscInt n,i,jrow,j,*ridx=PETSC_NULL; 1040 PetscScalar sum; 1041 PetscTruth usecprow=a->compressedrow.use; 1042 #endif 1043 1044 PetscFunctionBegin; 1045 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1046 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1047 if (zz != yy) { 1048 ierr = VecGetArray(zz,&z);CHKERRQ(ierr); 1049 } else { 1050 z = y; 1051 } 1052 1053 aj = a->j; 1054 aa = a->a; 1055 ii = a->i; 1056 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 1057 fortranmultaddaij_(&m,x,ii,aj,aa,y,z); 1058 #else 1059 if (usecprow){ /* use compressed row format */ 1060 if (zz != yy){ 1061 ierr = PetscMemcpy(z,y,m*sizeof(PetscScalar));CHKERRQ(ierr); 1062 } 1063 m = a->compressedrow.nrows; 1064 ii = a->compressedrow.i; 1065 ridx = a->compressedrow.rindex; 1066 for (i=0; i<m; i++){ 1067 n = ii[i+1] - ii[i]; 1068 aj = a->j + ii[i]; 1069 aa = a->a + ii[i]; 1070 sum = y[*ridx]; 1071 for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; 1072 z[*ridx++] = sum; 1073 } 1074 } else { /* do not use compressed row format */ 1075 for (i=0; i<m; i++) { 1076 jrow = ii[i]; 1077 n = ii[i+1] - jrow; 1078 sum = y[i]; 1079 for (j=0; j<n; j++) { 1080 sum += aa[jrow]*x[aj[jrow]]; jrow++; 1081 } 1082 z[i] = sum; 1083 } 1084 } 1085 #endif 1086 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1087 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1088 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1089 if (zz != yy) { 1090 ierr = VecRestoreArray(zz,&z);CHKERRQ(ierr); 1091 } 1092 PetscFunctionReturn(0); 1093 } 1094 1095 /* 1096 Adds diagonal pointers to sparse matrix structure. 1097 */ 1098 #undef __FUNCT__ 1099 #define __FUNCT__ "MatMarkDiagonal_SeqAIJ" 1100 PetscErrorCode MatMarkDiagonal_SeqAIJ(Mat A) 1101 { 1102 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1103 PetscErrorCode ierr; 1104 PetscInt i,j,m = A->rmap->n; 1105 1106 PetscFunctionBegin; 1107 if (!a->diag) { 1108 ierr = PetscMalloc(m*sizeof(PetscInt),&a->diag);CHKERRQ(ierr); 1109 ierr = PetscLogObjectMemory(A, m*sizeof(PetscInt));CHKERRQ(ierr); 1110 } 1111 for (i=0; i<A->rmap->n; i++) { 1112 a->diag[i] = a->i[i+1]; 1113 for (j=a->i[i]; j<a->i[i+1]; j++) { 1114 if (a->j[j] == i) { 1115 a->diag[i] = j; 1116 break; 1117 } 1118 } 1119 } 1120 PetscFunctionReturn(0); 1121 } 1122 1123 /* 1124 Checks for missing diagonals 1125 */ 1126 #undef __FUNCT__ 1127 #define __FUNCT__ "MatMissingDiagonal_SeqAIJ" 1128 PetscErrorCode MatMissingDiagonal_SeqAIJ(Mat A,PetscTruth *missing,PetscInt *d) 1129 { 1130 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1131 PetscInt *diag,*jj = a->j,i; 1132 1133 PetscFunctionBegin; 1134 *missing = PETSC_FALSE; 1135 if (A->rmap->n > 0 && !jj) { 1136 *missing = PETSC_TRUE; 1137 if (d) *d = 0; 1138 PetscInfo(A,"Matrix has no entries therefor is missing diagonal"); 1139 } else { 1140 diag = a->diag; 1141 for (i=0; i<A->rmap->n; i++) { 1142 if (jj[diag[i]] != i) { 1143 *missing = PETSC_TRUE; 1144 if (d) *d = i; 1145 PetscInfo1(A,"Matrix is missing diagonal number %D",i); 1146 } 1147 } 1148 } 1149 PetscFunctionReturn(0); 1150 } 1151 1152 EXTERN_C_BEGIN 1153 #undef __FUNCT__ 1154 #define __FUNCT__ "MatInvertDiagonal_SeqAIJ" 1155 PetscErrorCode PETSCMAT_DLLEXPORT MatInvertDiagonal_SeqAIJ(Mat A,PetscScalar omega,PetscScalar fshift) 1156 { 1157 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 1158 PetscErrorCode ierr; 1159 PetscInt i,*diag,m = A->rmap->n; 1160 MatScalar *v = a->a; 1161 PetscScalar *idiag,*mdiag; 1162 1163 PetscFunctionBegin; 1164 if (a->idiagvalid) PetscFunctionReturn(0); 1165 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1166 diag = a->diag; 1167 if (!a->idiag) { 1168 ierr = PetscMalloc3(m,PetscScalar,&a->idiag,m,PetscScalar,&a->mdiag,m,PetscScalar,&a->ssor_work);CHKERRQ(ierr); 1169 ierr = PetscLogObjectMemory(A, 3*m*sizeof(PetscScalar));CHKERRQ(ierr); 1170 v = a->a; 1171 } 1172 mdiag = a->mdiag; 1173 idiag = a->idiag; 1174 1175 if (omega == 1.0 && !PetscAbsScalar(fshift)) { 1176 for (i=0; i<m; i++) { 1177 mdiag[i] = v[diag[i]]; 1178 if (!PetscAbsScalar(mdiag[i])) SETERRQ1(PETSC_ERR_ARG_INCOMP,"Zero diagonal on row %D",i); 1179 idiag[i] = 1.0/v[diag[i]]; 1180 } 1181 ierr = PetscLogFlops(m);CHKERRQ(ierr); 1182 } else { 1183 for (i=0; i<m; i++) { 1184 mdiag[i] = v[diag[i]]; 1185 idiag[i] = omega/(fshift + v[diag[i]]); 1186 } 1187 ierr = PetscLogFlops(2.0*m);CHKERRQ(ierr); 1188 } 1189 a->idiagvalid = PETSC_TRUE; 1190 PetscFunctionReturn(0); 1191 } 1192 EXTERN_C_END 1193 1194 #include "../src/mat/impls/aij/seq/ftn-kernels/frelax.h" 1195 #undef __FUNCT__ 1196 #define __FUNCT__ "MatSOR_SeqAIJ" 1197 PetscErrorCode MatSOR_SeqAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx) 1198 { 1199 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1200 PetscScalar *x,d,sum,*t,scale; 1201 const MatScalar *v = a->a,*idiag=0,*mdiag; 1202 const PetscScalar *b, *bs,*xb, *ts; 1203 PetscErrorCode ierr; 1204 PetscInt n = A->cmap->n,m = A->rmap->n,i; 1205 const PetscInt *idx,*diag; 1206 1207 PetscFunctionBegin; 1208 its = its*lits; 1209 1210 if (fshift != a->fshift || omega != a->omega) a->idiagvalid = PETSC_FALSE; /* must recompute idiag[] */ 1211 if (!a->idiagvalid) {ierr = MatInvertDiagonal_SeqAIJ(A,omega,fshift);CHKERRQ(ierr);} 1212 a->fshift = fshift; 1213 a->omega = omega; 1214 1215 diag = a->diag; 1216 t = a->ssor_work; 1217 idiag = a->idiag; 1218 mdiag = a->mdiag; 1219 1220 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1221 if (xx != bb) { 1222 ierr = VecGetArray(bb,(PetscScalar**)&b);CHKERRQ(ierr); 1223 } else { 1224 b = x; 1225 } 1226 CHKMEMQ; 1227 /* We count flops by assuming the upper triangular and lower triangular parts have the same number of nonzeros */ 1228 if (flag == SOR_APPLY_UPPER) { 1229 /* apply (U + D/omega) to the vector */ 1230 bs = b; 1231 for (i=0; i<m; i++) { 1232 d = fshift + mdiag[i]; 1233 n = a->i[i+1] - diag[i] - 1; 1234 idx = a->j + diag[i] + 1; 1235 v = a->a + diag[i] + 1; 1236 sum = b[i]*d/omega; 1237 PetscSparseDensePlusDot(sum,bs,v,idx,n); 1238 x[i] = sum; 1239 } 1240 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1241 if (bb != xx) {ierr = VecRestoreArray(bb,(PetscScalar**)&b);CHKERRQ(ierr);} 1242 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1243 PetscFunctionReturn(0); 1244 } 1245 1246 if (flag == SOR_APPLY_LOWER) { 1247 SETERRQ(PETSC_ERR_SUP,"SOR_APPLY_LOWER is not implemented"); 1248 } else if (flag & SOR_EISENSTAT) { 1249 /* Let A = L + U + D; where L is lower trianglar, 1250 U is upper triangular, E = D/omega; This routine applies 1251 1252 (L + E)^{-1} A (U + E)^{-1} 1253 1254 to a vector efficiently using Eisenstat's trick. 1255 */ 1256 scale = (2.0/omega) - 1.0; 1257 1258 /* x = (E + U)^{-1} b */ 1259 for (i=m-1; i>=0; i--) { 1260 n = a->i[i+1] - diag[i] - 1; 1261 idx = a->j + diag[i] + 1; 1262 v = a->a + diag[i] + 1; 1263 sum = b[i]; 1264 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1265 x[i] = sum*idiag[i]; 1266 } 1267 1268 /* t = b - (2*E - D)x */ 1269 v = a->a; 1270 for (i=0; i<m; i++) { t[i] = b[i] - scale*(v[*diag++])*x[i]; } 1271 1272 /* t = (E + L)^{-1}t */ 1273 ts = t; 1274 diag = a->diag; 1275 for (i=0; i<m; i++) { 1276 n = diag[i] - a->i[i]; 1277 idx = a->j + a->i[i]; 1278 v = a->a + a->i[i]; 1279 sum = t[i]; 1280 PetscSparseDenseMinusDot(sum,ts,v,idx,n); 1281 t[i] = sum*idiag[i]; 1282 /* x = x + t */ 1283 x[i] += t[i]; 1284 } 1285 1286 ierr = PetscLogFlops(6.0*m-1 + 2.0*a->nz);CHKERRQ(ierr); 1287 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1288 if (bb != xx) {ierr = VecRestoreArray(bb,(PetscScalar**)&b);CHKERRQ(ierr);} 1289 PetscFunctionReturn(0); 1290 } 1291 if (flag & SOR_ZERO_INITIAL_GUESS) { 1292 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1293 for (i=0; i<m; i++) { 1294 n = diag[i] - a->i[i]; 1295 idx = a->j + a->i[i]; 1296 v = a->a + a->i[i]; 1297 sum = b[i]; 1298 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1299 t[i] = sum; 1300 x[i] = sum*idiag[i]; 1301 } 1302 xb = t; 1303 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1304 } else xb = b; 1305 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1306 for (i=m-1; i>=0; i--) { 1307 n = a->i[i+1] - diag[i] - 1; 1308 idx = a->j + diag[i] + 1; 1309 v = a->a + diag[i] + 1; 1310 sum = xb[i]; 1311 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1312 if (xb == b) { 1313 x[i] = sum*idiag[i]; 1314 } else { 1315 x[i] = (1-omega)*x[i] + sum*idiag[i]; 1316 } 1317 } 1318 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1319 } 1320 its--; 1321 } 1322 while (its--) { 1323 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 1324 for (i=0; i<m; i++) { 1325 n = a->i[i+1] - a->i[i]; 1326 idx = a->j + a->i[i]; 1327 v = a->a + a->i[i]; 1328 sum = b[i]; 1329 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1330 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1331 } 1332 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1333 } 1334 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 1335 for (i=m-1; i>=0; i--) { 1336 n = a->i[i+1] - a->i[i]; 1337 idx = a->j + a->i[i]; 1338 v = a->a + a->i[i]; 1339 sum = b[i]; 1340 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1341 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1342 } 1343 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1344 } 1345 } 1346 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1347 if (bb != xx) {ierr = VecRestoreArray(bb,(PetscScalar**)&b);CHKERRQ(ierr);} 1348 CHKMEMQ; PetscFunctionReturn(0); 1349 } 1350 1351 1352 #undef __FUNCT__ 1353 #define __FUNCT__ "MatGetInfo_SeqAIJ" 1354 PetscErrorCode MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1355 { 1356 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1357 1358 PetscFunctionBegin; 1359 info->block_size = 1.0; 1360 info->nz_allocated = (double)a->maxnz; 1361 info->nz_used = (double)a->nz; 1362 info->nz_unneeded = (double)(a->maxnz - a->nz); 1363 info->assemblies = (double)A->num_ass; 1364 info->mallocs = (double)a->reallocs; 1365 info->memory = ((PetscObject)A)->mem; 1366 if (A->factor) { 1367 info->fill_ratio_given = A->info.fill_ratio_given; 1368 info->fill_ratio_needed = A->info.fill_ratio_needed; 1369 info->factor_mallocs = A->info.factor_mallocs; 1370 } else { 1371 info->fill_ratio_given = 0; 1372 info->fill_ratio_needed = 0; 1373 info->factor_mallocs = 0; 1374 } 1375 PetscFunctionReturn(0); 1376 } 1377 1378 #undef __FUNCT__ 1379 #define __FUNCT__ "MatZeroRows_SeqAIJ" 1380 PetscErrorCode MatZeroRows_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag) 1381 { 1382 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1383 PetscInt i,m = A->rmap->n - 1,d; 1384 PetscErrorCode ierr; 1385 PetscTruth missing; 1386 1387 PetscFunctionBegin; 1388 if (a->keepnonzeropattern) { 1389 for (i=0; i<N; i++) { 1390 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1391 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 1392 } 1393 if (diag != 0.0) { 1394 ierr = MatMissingDiagonal_SeqAIJ(A,&missing,&d);CHKERRQ(ierr); 1395 if (missing) SETERRQ1(PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in row %D",d); 1396 for (i=0; i<N; i++) { 1397 a->a[a->diag[rows[i]]] = diag; 1398 } 1399 } 1400 A->same_nonzero = PETSC_TRUE; 1401 } else { 1402 if (diag != 0.0) { 1403 for (i=0; i<N; i++) { 1404 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1405 if (a->ilen[rows[i]] > 0) { 1406 a->ilen[rows[i]] = 1; 1407 a->a[a->i[rows[i]]] = diag; 1408 a->j[a->i[rows[i]]] = rows[i]; 1409 } else { /* in case row was completely empty */ 1410 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES);CHKERRQ(ierr); 1411 } 1412 } 1413 } else { 1414 for (i=0; i<N; i++) { 1415 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1416 a->ilen[rows[i]] = 0; 1417 } 1418 } 1419 A->same_nonzero = PETSC_FALSE; 1420 } 1421 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1422 PetscFunctionReturn(0); 1423 } 1424 1425 #undef __FUNCT__ 1426 #define __FUNCT__ "MatGetRow_SeqAIJ" 1427 PetscErrorCode MatGetRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1428 { 1429 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1430 PetscInt *itmp; 1431 1432 PetscFunctionBegin; 1433 if (row < 0 || row >= A->rmap->n) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Row %D out of range",row); 1434 1435 *nz = a->i[row+1] - a->i[row]; 1436 if (v) *v = a->a + a->i[row]; 1437 if (idx) { 1438 itmp = a->j + a->i[row]; 1439 if (*nz) { 1440 *idx = itmp; 1441 } 1442 else *idx = 0; 1443 } 1444 PetscFunctionReturn(0); 1445 } 1446 1447 /* remove this function? */ 1448 #undef __FUNCT__ 1449 #define __FUNCT__ "MatRestoreRow_SeqAIJ" 1450 PetscErrorCode MatRestoreRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1451 { 1452 PetscFunctionBegin; 1453 PetscFunctionReturn(0); 1454 } 1455 1456 #undef __FUNCT__ 1457 #define __FUNCT__ "MatNorm_SeqAIJ" 1458 PetscErrorCode MatNorm_SeqAIJ(Mat A,NormType type,PetscReal *nrm) 1459 { 1460 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1461 MatScalar *v = a->a; 1462 PetscReal sum = 0.0; 1463 PetscErrorCode ierr; 1464 PetscInt i,j; 1465 1466 PetscFunctionBegin; 1467 if (type == NORM_FROBENIUS) { 1468 for (i=0; i<a->nz; i++) { 1469 #if defined(PETSC_USE_COMPLEX) 1470 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1471 #else 1472 sum += (*v)*(*v); v++; 1473 #endif 1474 } 1475 *nrm = sqrt(sum); 1476 } else if (type == NORM_1) { 1477 PetscReal *tmp; 1478 PetscInt *jj = a->j; 1479 ierr = PetscMalloc((A->cmap->n+1)*sizeof(PetscReal),&tmp);CHKERRQ(ierr); 1480 ierr = PetscMemzero(tmp,A->cmap->n*sizeof(PetscReal));CHKERRQ(ierr); 1481 *nrm = 0.0; 1482 for (j=0; j<a->nz; j++) { 1483 tmp[*jj++] += PetscAbsScalar(*v); v++; 1484 } 1485 for (j=0; j<A->cmap->n; j++) { 1486 if (tmp[j] > *nrm) *nrm = tmp[j]; 1487 } 1488 ierr = PetscFree(tmp);CHKERRQ(ierr); 1489 } else if (type == NORM_INFINITY) { 1490 *nrm = 0.0; 1491 for (j=0; j<A->rmap->n; j++) { 1492 v = a->a + a->i[j]; 1493 sum = 0.0; 1494 for (i=0; i<a->i[j+1]-a->i[j]; i++) { 1495 sum += PetscAbsScalar(*v); v++; 1496 } 1497 if (sum > *nrm) *nrm = sum; 1498 } 1499 } else { 1500 SETERRQ(PETSC_ERR_SUP,"No support for two norm"); 1501 } 1502 PetscFunctionReturn(0); 1503 } 1504 1505 #undef __FUNCT__ 1506 #define __FUNCT__ "MatTranspose_SeqAIJ" 1507 PetscErrorCode MatTranspose_SeqAIJ(Mat A,MatReuse reuse,Mat *B) 1508 { 1509 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1510 Mat C; 1511 PetscErrorCode ierr; 1512 PetscInt i,*aj = a->j,*ai = a->i,m = A->rmap->n,len,*col; 1513 MatScalar *array = a->a; 1514 1515 PetscFunctionBegin; 1516 if (reuse == MAT_REUSE_MATRIX && A == *B && m != A->cmap->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Square matrix only for in-place"); 1517 1518 if (reuse == MAT_INITIAL_MATRIX || *B == A) { 1519 ierr = PetscMalloc((1+A->cmap->n)*sizeof(PetscInt),&col);CHKERRQ(ierr); 1520 ierr = PetscMemzero(col,(1+A->cmap->n)*sizeof(PetscInt));CHKERRQ(ierr); 1521 1522 for (i=0; i<ai[m]; i++) col[aj[i]] += 1; 1523 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 1524 ierr = MatSetSizes(C,A->cmap->n,m,A->cmap->n,m);CHKERRQ(ierr); 1525 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 1526 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,col);CHKERRQ(ierr); 1527 ierr = PetscFree(col);CHKERRQ(ierr); 1528 } else { 1529 C = *B; 1530 } 1531 1532 for (i=0; i<m; i++) { 1533 len = ai[i+1]-ai[i]; 1534 ierr = MatSetValues_SeqAIJ(C,len,aj,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 1535 array += len; 1536 aj += len; 1537 } 1538 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1539 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1540 1541 if (reuse == MAT_INITIAL_MATRIX || *B != A) { 1542 *B = C; 1543 } else { 1544 ierr = MatHeaderCopy(A,C);CHKERRQ(ierr); 1545 } 1546 PetscFunctionReturn(0); 1547 } 1548 1549 EXTERN_C_BEGIN 1550 #undef __FUNCT__ 1551 #define __FUNCT__ "MatIsTranspose_SeqAIJ" 1552 PetscErrorCode PETSCMAT_DLLEXPORT MatIsTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscTruth *f) 1553 { 1554 Mat_SeqAIJ *aij = (Mat_SeqAIJ *) A->data,*bij = (Mat_SeqAIJ*) A->data; 1555 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 1556 MatScalar *va,*vb; 1557 PetscErrorCode ierr; 1558 PetscInt ma,na,mb,nb, i; 1559 1560 PetscFunctionBegin; 1561 bij = (Mat_SeqAIJ *) B->data; 1562 1563 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 1564 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 1565 if (ma!=nb || na!=mb){ 1566 *f = PETSC_FALSE; 1567 PetscFunctionReturn(0); 1568 } 1569 aii = aij->i; bii = bij->i; 1570 adx = aij->j; bdx = bij->j; 1571 va = aij->a; vb = bij->a; 1572 ierr = PetscMalloc(ma*sizeof(PetscInt),&aptr);CHKERRQ(ierr); 1573 ierr = PetscMalloc(mb*sizeof(PetscInt),&bptr);CHKERRQ(ierr); 1574 for (i=0; i<ma; i++) aptr[i] = aii[i]; 1575 for (i=0; i<mb; i++) bptr[i] = bii[i]; 1576 1577 *f = PETSC_TRUE; 1578 for (i=0; i<ma; i++) { 1579 while (aptr[i]<aii[i+1]) { 1580 PetscInt idc,idr; 1581 PetscScalar vc,vr; 1582 /* column/row index/value */ 1583 idc = adx[aptr[i]]; 1584 idr = bdx[bptr[idc]]; 1585 vc = va[aptr[i]]; 1586 vr = vb[bptr[idc]]; 1587 if (i!=idr || PetscAbsScalar(vc-vr) > tol) { 1588 *f = PETSC_FALSE; 1589 goto done; 1590 } else { 1591 aptr[i]++; 1592 if (B || i!=idc) bptr[idc]++; 1593 } 1594 } 1595 } 1596 done: 1597 ierr = PetscFree(aptr);CHKERRQ(ierr); 1598 if (B) { 1599 ierr = PetscFree(bptr);CHKERRQ(ierr); 1600 } 1601 PetscFunctionReturn(0); 1602 } 1603 EXTERN_C_END 1604 1605 EXTERN_C_BEGIN 1606 #undef __FUNCT__ 1607 #define __FUNCT__ "MatIsHermitianTranspose_SeqAIJ" 1608 PetscErrorCode PETSCMAT_DLLEXPORT MatIsHermitianTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscTruth *f) 1609 { 1610 Mat_SeqAIJ *aij = (Mat_SeqAIJ *) A->data,*bij = (Mat_SeqAIJ*) A->data; 1611 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 1612 MatScalar *va,*vb; 1613 PetscErrorCode ierr; 1614 PetscInt ma,na,mb,nb, i; 1615 1616 PetscFunctionBegin; 1617 bij = (Mat_SeqAIJ *) B->data; 1618 1619 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 1620 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 1621 if (ma!=nb || na!=mb){ 1622 *f = PETSC_FALSE; 1623 PetscFunctionReturn(0); 1624 } 1625 aii = aij->i; bii = bij->i; 1626 adx = aij->j; bdx = bij->j; 1627 va = aij->a; vb = bij->a; 1628 ierr = PetscMalloc(ma*sizeof(PetscInt),&aptr);CHKERRQ(ierr); 1629 ierr = PetscMalloc(mb*sizeof(PetscInt),&bptr);CHKERRQ(ierr); 1630 for (i=0; i<ma; i++) aptr[i] = aii[i]; 1631 for (i=0; i<mb; i++) bptr[i] = bii[i]; 1632 1633 *f = PETSC_TRUE; 1634 for (i=0; i<ma; i++) { 1635 while (aptr[i]<aii[i+1]) { 1636 PetscInt idc,idr; 1637 PetscScalar vc,vr; 1638 /* column/row index/value */ 1639 idc = adx[aptr[i]]; 1640 idr = bdx[bptr[idc]]; 1641 vc = va[aptr[i]]; 1642 vr = vb[bptr[idc]]; 1643 if (i!=idr || PetscAbsScalar(vc-PetscConj(vr)) > tol) { 1644 *f = PETSC_FALSE; 1645 goto done; 1646 } else { 1647 aptr[i]++; 1648 if (B || i!=idc) bptr[idc]++; 1649 } 1650 } 1651 } 1652 done: 1653 ierr = PetscFree(aptr);CHKERRQ(ierr); 1654 if (B) { 1655 ierr = PetscFree(bptr);CHKERRQ(ierr); 1656 } 1657 PetscFunctionReturn(0); 1658 } 1659 EXTERN_C_END 1660 1661 #undef __FUNCT__ 1662 #define __FUNCT__ "MatIsSymmetric_SeqAIJ" 1663 PetscErrorCode MatIsSymmetric_SeqAIJ(Mat A,PetscReal tol,PetscTruth *f) 1664 { 1665 PetscErrorCode ierr; 1666 PetscFunctionBegin; 1667 ierr = MatIsTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 1668 PetscFunctionReturn(0); 1669 } 1670 1671 #undef __FUNCT__ 1672 #define __FUNCT__ "MatIsHermitian_SeqAIJ" 1673 PetscErrorCode MatIsHermitian_SeqAIJ(Mat A,PetscReal tol,PetscTruth *f) 1674 { 1675 PetscErrorCode ierr; 1676 PetscFunctionBegin; 1677 ierr = MatIsHermitianTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 1678 PetscFunctionReturn(0); 1679 } 1680 1681 #undef __FUNCT__ 1682 #define __FUNCT__ "MatDiagonalScale_SeqAIJ" 1683 PetscErrorCode MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 1684 { 1685 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1686 PetscScalar *l,*r,x; 1687 MatScalar *v; 1688 PetscErrorCode ierr; 1689 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,M,nz = a->nz,*jj; 1690 1691 PetscFunctionBegin; 1692 if (ll) { 1693 /* The local size is used so that VecMPI can be passed to this routine 1694 by MatDiagonalScale_MPIAIJ */ 1695 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 1696 if (m != A->rmap->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 1697 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 1698 v = a->a; 1699 for (i=0; i<m; i++) { 1700 x = l[i]; 1701 M = a->i[i+1] - a->i[i]; 1702 for (j=0; j<M; j++) { (*v++) *= x;} 1703 } 1704 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 1705 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 1706 } 1707 if (rr) { 1708 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 1709 if (n != A->cmap->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 1710 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 1711 v = a->a; jj = a->j; 1712 for (i=0; i<nz; i++) { 1713 (*v++) *= r[*jj++]; 1714 } 1715 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 1716 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 1717 } 1718 PetscFunctionReturn(0); 1719 } 1720 1721 #undef __FUNCT__ 1722 #define __FUNCT__ "MatGetSubMatrix_SeqAIJ" 1723 PetscErrorCode MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,PetscInt csize,MatReuse scall,Mat *B) 1724 { 1725 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 1726 PetscErrorCode ierr; 1727 PetscInt *smap,i,k,kstart,kend,oldcols = A->cmap->n,*lens; 1728 PetscInt row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 1729 const PetscInt *irow,*icol; 1730 PetscInt nrows,ncols; 1731 PetscInt *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 1732 MatScalar *a_new,*mat_a; 1733 Mat C; 1734 PetscTruth stride,sorted; 1735 1736 PetscFunctionBegin; 1737 ierr = ISSorted(isrow,&sorted);CHKERRQ(ierr); 1738 if (!sorted) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"ISrow is not sorted"); 1739 ierr = ISSorted(iscol,&sorted);CHKERRQ(ierr); 1740 if (!sorted) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"IScol is not sorted"); 1741 1742 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 1743 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 1744 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 1745 1746 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 1747 ierr = ISStride(iscol,&stride);CHKERRQ(ierr); 1748 if (stride && step == 1) { 1749 /* special case of contiguous rows */ 1750 ierr = PetscMalloc2(nrows,PetscInt,&lens,nrows,PetscInt,&starts);CHKERRQ(ierr); 1751 /* loop over new rows determining lens and starting points */ 1752 for (i=0; i<nrows; i++) { 1753 kstart = ai[irow[i]]; 1754 kend = kstart + ailen[irow[i]]; 1755 for (k=kstart; k<kend; k++) { 1756 if (aj[k] >= first) { 1757 starts[i] = k; 1758 break; 1759 } 1760 } 1761 sum = 0; 1762 while (k < kend) { 1763 if (aj[k++] >= first+ncols) break; 1764 sum++; 1765 } 1766 lens[i] = sum; 1767 } 1768 /* create submatrix */ 1769 if (scall == MAT_REUSE_MATRIX) { 1770 PetscInt n_cols,n_rows; 1771 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1772 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 1773 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 1774 C = *B; 1775 } else { 1776 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 1777 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1778 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 1779 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 1780 } 1781 c = (Mat_SeqAIJ*)C->data; 1782 1783 /* loop over rows inserting into submatrix */ 1784 a_new = c->a; 1785 j_new = c->j; 1786 i_new = c->i; 1787 1788 for (i=0; i<nrows; i++) { 1789 ii = starts[i]; 1790 lensi = lens[i]; 1791 for (k=0; k<lensi; k++) { 1792 *j_new++ = aj[ii+k] - first; 1793 } 1794 ierr = PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(PetscScalar));CHKERRQ(ierr); 1795 a_new += lensi; 1796 i_new[i+1] = i_new[i] + lensi; 1797 c->ilen[i] = lensi; 1798 } 1799 ierr = PetscFree2(lens,starts);CHKERRQ(ierr); 1800 } else { 1801 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1802 ierr = PetscMalloc(oldcols*sizeof(PetscInt),&smap);CHKERRQ(ierr); 1803 ierr = PetscMemzero(smap,oldcols*sizeof(PetscInt));CHKERRQ(ierr); 1804 ierr = PetscMalloc((1+nrows)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 1805 for (i=0; i<ncols; i++) smap[icol[i]] = i+1; 1806 /* determine lens of each row */ 1807 for (i=0; i<nrows; i++) { 1808 kstart = ai[irow[i]]; 1809 kend = kstart + a->ilen[irow[i]]; 1810 lens[i] = 0; 1811 for (k=kstart; k<kend; k++) { 1812 if (smap[aj[k]]) { 1813 lens[i]++; 1814 } 1815 } 1816 } 1817 /* Create and fill new matrix */ 1818 if (scall == MAT_REUSE_MATRIX) { 1819 PetscTruth equal; 1820 1821 c = (Mat_SeqAIJ *)((*B)->data); 1822 if ((*B)->rmap->n != nrows || (*B)->cmap->n != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 1823 ierr = PetscMemcmp(c->ilen,lens,(*B)->rmap->n*sizeof(PetscInt),&equal);CHKERRQ(ierr); 1824 if (!equal) { 1825 SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 1826 } 1827 ierr = PetscMemzero(c->ilen,(*B)->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 1828 C = *B; 1829 } else { 1830 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 1831 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1832 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 1833 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 1834 } 1835 c = (Mat_SeqAIJ *)(C->data); 1836 for (i=0; i<nrows; i++) { 1837 row = irow[i]; 1838 kstart = ai[row]; 1839 kend = kstart + a->ilen[row]; 1840 mat_i = c->i[i]; 1841 mat_j = c->j + mat_i; 1842 mat_a = c->a + mat_i; 1843 mat_ilen = c->ilen + i; 1844 for (k=kstart; k<kend; k++) { 1845 if ((tcol=smap[a->j[k]])) { 1846 *mat_j++ = tcol - 1; 1847 *mat_a++ = a->a[k]; 1848 (*mat_ilen)++; 1849 1850 } 1851 } 1852 } 1853 /* Free work space */ 1854 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1855 ierr = PetscFree(smap);CHKERRQ(ierr); 1856 ierr = PetscFree(lens);CHKERRQ(ierr); 1857 } 1858 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1859 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1860 1861 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1862 *B = C; 1863 PetscFunctionReturn(0); 1864 } 1865 1866 #undef __FUNCT__ 1867 #define __FUNCT__ "MatILUFactor_SeqAIJ" 1868 PetscErrorCode MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,const MatFactorInfo *info) 1869 { 1870 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1871 PetscErrorCode ierr; 1872 Mat outA; 1873 PetscTruth row_identity,col_identity; 1874 1875 PetscFunctionBegin; 1876 if (info->levels != 0) SETERRQ(PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 1877 1878 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 1879 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 1880 1881 outA = inA; 1882 outA->factor = MAT_FACTOR_LU; 1883 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1884 if (a->row) { ierr = ISDestroy(a->row);CHKERRQ(ierr);} 1885 a->row = row; 1886 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 1887 if (a->col) { ierr = ISDestroy(a->col);CHKERRQ(ierr);} 1888 a->col = col; 1889 1890 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 1891 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} /* need to remove old one */ 1892 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 1893 ierr = PetscLogObjectParent(inA,a->icol);CHKERRQ(ierr); 1894 1895 if (!a->solve_work) { /* this matrix may have been factored before */ 1896 ierr = PetscMalloc((inA->rmap->n+1)*sizeof(PetscScalar),&a->solve_work);CHKERRQ(ierr); 1897 ierr = PetscLogObjectMemory(inA, (inA->rmap->n+1)*sizeof(PetscScalar));CHKERRQ(ierr); 1898 } 1899 1900 ierr = MatMarkDiagonal_SeqAIJ(inA);CHKERRQ(ierr); 1901 if (row_identity && col_identity) { 1902 ierr = MatLUFactorNumeric_SeqAIJ_inplace(outA,inA,info);CHKERRQ(ierr); 1903 } else { 1904 ierr = MatLUFactorNumeric_SeqAIJ_InplaceWithPerm(outA,inA,info);CHKERRQ(ierr); 1905 } 1906 PetscFunctionReturn(0); 1907 } 1908 1909 #include "petscblaslapack.h" 1910 #undef __FUNCT__ 1911 #define __FUNCT__ "MatScale_SeqAIJ" 1912 PetscErrorCode MatScale_SeqAIJ(Mat inA,PetscScalar alpha) 1913 { 1914 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1915 PetscScalar oalpha = alpha; 1916 PetscErrorCode ierr; 1917 PetscBLASInt one = 1,bnz = PetscBLASIntCast(a->nz); 1918 1919 PetscFunctionBegin; 1920 BLASscal_(&bnz,&oalpha,a->a,&one); 1921 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1922 PetscFunctionReturn(0); 1923 } 1924 1925 #undef __FUNCT__ 1926 #define __FUNCT__ "MatGetSubMatrices_SeqAIJ" 1927 PetscErrorCode MatGetSubMatrices_SeqAIJ(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1928 { 1929 PetscErrorCode ierr; 1930 PetscInt i; 1931 1932 PetscFunctionBegin; 1933 if (scall == MAT_INITIAL_MATRIX) { 1934 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1935 } 1936 1937 for (i=0; i<n; i++) { 1938 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1939 } 1940 PetscFunctionReturn(0); 1941 } 1942 1943 #undef __FUNCT__ 1944 #define __FUNCT__ "MatIncreaseOverlap_SeqAIJ" 1945 PetscErrorCode MatIncreaseOverlap_SeqAIJ(Mat A,PetscInt is_max,IS is[],PetscInt ov) 1946 { 1947 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1948 PetscErrorCode ierr; 1949 PetscInt row,i,j,k,l,m,n,*nidx,isz,val; 1950 const PetscInt *idx; 1951 PetscInt start,end,*ai,*aj; 1952 PetscBT table; 1953 1954 PetscFunctionBegin; 1955 m = A->rmap->n; 1956 ai = a->i; 1957 aj = a->j; 1958 1959 if (ov < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"illegal negative overlap value used"); 1960 1961 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&nidx);CHKERRQ(ierr); 1962 ierr = PetscBTCreate(m,table);CHKERRQ(ierr); 1963 1964 for (i=0; i<is_max; i++) { 1965 /* Initialize the two local arrays */ 1966 isz = 0; 1967 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 1968 1969 /* Extract the indices, assume there can be duplicate entries */ 1970 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 1971 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 1972 1973 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 1974 for (j=0; j<n ; ++j){ 1975 if(!PetscBTLookupSet(table,idx[j])) { nidx[isz++] = idx[j];} 1976 } 1977 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 1978 ierr = ISDestroy(is[i]);CHKERRQ(ierr); 1979 1980 k = 0; 1981 for (j=0; j<ov; j++){ /* for each overlap */ 1982 n = isz; 1983 for (; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 1984 row = nidx[k]; 1985 start = ai[row]; 1986 end = ai[row+1]; 1987 for (l = start; l<end ; l++){ 1988 val = aj[l] ; 1989 if (!PetscBTLookupSet(table,val)) {nidx[isz++] = val;} 1990 } 1991 } 1992 } 1993 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,(is+i));CHKERRQ(ierr); 1994 } 1995 ierr = PetscBTDestroy(table);CHKERRQ(ierr); 1996 ierr = PetscFree(nidx);CHKERRQ(ierr); 1997 PetscFunctionReturn(0); 1998 } 1999 2000 /* -------------------------------------------------------------- */ 2001 #undef __FUNCT__ 2002 #define __FUNCT__ "MatPermute_SeqAIJ" 2003 PetscErrorCode MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 2004 { 2005 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2006 PetscErrorCode ierr; 2007 PetscInt i,nz,m = A->rmap->n,n = A->cmap->n; 2008 const PetscInt *row,*col; 2009 PetscInt *cnew,j,*lens; 2010 IS icolp,irowp; 2011 PetscInt *cwork; 2012 PetscScalar *vwork; 2013 2014 PetscFunctionBegin; 2015 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 2016 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 2017 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 2018 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 2019 2020 /* determine lengths of permuted rows */ 2021 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 2022 for (i=0; i<m; i++) { 2023 lens[row[i]] = a->i[i+1] - a->i[i]; 2024 } 2025 ierr = MatCreate(((PetscObject)A)->comm,B);CHKERRQ(ierr); 2026 ierr = MatSetSizes(*B,m,n,m,n);CHKERRQ(ierr); 2027 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 2028 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*B,0,lens);CHKERRQ(ierr); 2029 ierr = PetscFree(lens);CHKERRQ(ierr); 2030 2031 ierr = PetscMalloc(n*sizeof(PetscInt),&cnew);CHKERRQ(ierr); 2032 for (i=0; i<m; i++) { 2033 ierr = MatGetRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2034 for (j=0; j<nz; j++) { cnew[j] = col[cwork[j]];} 2035 ierr = MatSetValues_SeqAIJ(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 2036 ierr = MatRestoreRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2037 } 2038 ierr = PetscFree(cnew);CHKERRQ(ierr); 2039 (*B)->assembled = PETSC_FALSE; 2040 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2041 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2042 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 2043 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 2044 ierr = ISDestroy(irowp);CHKERRQ(ierr); 2045 ierr = ISDestroy(icolp);CHKERRQ(ierr); 2046 PetscFunctionReturn(0); 2047 } 2048 2049 #undef __FUNCT__ 2050 #define __FUNCT__ "MatCopy_SeqAIJ" 2051 PetscErrorCode MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 2052 { 2053 PetscErrorCode ierr; 2054 2055 PetscFunctionBegin; 2056 /* If the two matrices have the same copy implementation, use fast copy. */ 2057 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 2058 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2059 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 2060 2061 if (a->i[A->rmap->n] != b->i[B->rmap->n]) { 2062 SETERRQ(PETSC_ERR_ARG_INCOMP,"Number of nonzeros in two matrices are different"); 2063 } 2064 ierr = PetscMemcpy(b->a,a->a,(a->i[A->rmap->n])*sizeof(PetscScalar));CHKERRQ(ierr); 2065 } else { 2066 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 2067 } 2068 PetscFunctionReturn(0); 2069 } 2070 2071 #undef __FUNCT__ 2072 #define __FUNCT__ "MatSetUpPreallocation_SeqAIJ" 2073 PetscErrorCode MatSetUpPreallocation_SeqAIJ(Mat A) 2074 { 2075 PetscErrorCode ierr; 2076 2077 PetscFunctionBegin; 2078 ierr = MatSeqAIJSetPreallocation_SeqAIJ(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 2079 PetscFunctionReturn(0); 2080 } 2081 2082 #undef __FUNCT__ 2083 #define __FUNCT__ "MatGetArray_SeqAIJ" 2084 PetscErrorCode MatGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 2085 { 2086 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2087 PetscFunctionBegin; 2088 *array = a->a; 2089 PetscFunctionReturn(0); 2090 } 2091 2092 #undef __FUNCT__ 2093 #define __FUNCT__ "MatRestoreArray_SeqAIJ" 2094 PetscErrorCode MatRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 2095 { 2096 PetscFunctionBegin; 2097 PetscFunctionReturn(0); 2098 } 2099 2100 #undef __FUNCT__ 2101 #define __FUNCT__ "MatFDColoringApply_SeqAIJ" 2102 PetscErrorCode MatFDColoringApply_SeqAIJ(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx) 2103 { 2104 PetscErrorCode (*f)(void*,Vec,Vec,void*) = (PetscErrorCode (*)(void*,Vec,Vec,void *))coloring->f; 2105 PetscErrorCode ierr; 2106 PetscInt k,N,start,end,l,row,col,srow,**vscaleforrow,m1,m2; 2107 PetscScalar dx,*y,*xx,*w3_array; 2108 PetscScalar *vscale_array; 2109 PetscReal epsilon = coloring->error_rel,umin = coloring->umin; 2110 Vec w1,w2,w3; 2111 void *fctx = coloring->fctx; 2112 PetscTruth flg = PETSC_FALSE; 2113 2114 PetscFunctionBegin; 2115 if (!coloring->w1) { 2116 ierr = VecDuplicate(x1,&coloring->w1);CHKERRQ(ierr); 2117 ierr = PetscLogObjectParent(coloring,coloring->w1);CHKERRQ(ierr); 2118 ierr = VecDuplicate(x1,&coloring->w2);CHKERRQ(ierr); 2119 ierr = PetscLogObjectParent(coloring,coloring->w2);CHKERRQ(ierr); 2120 ierr = VecDuplicate(x1,&coloring->w3);CHKERRQ(ierr); 2121 ierr = PetscLogObjectParent(coloring,coloring->w3);CHKERRQ(ierr); 2122 } 2123 w1 = coloring->w1; w2 = coloring->w2; w3 = coloring->w3; 2124 2125 ierr = MatSetUnfactored(J);CHKERRQ(ierr); 2126 ierr = PetscOptionsGetTruth(((PetscObject)coloring)->prefix,"-mat_fd_coloring_dont_rezero",&flg,PETSC_NULL);CHKERRQ(ierr); 2127 if (flg) { 2128 ierr = PetscInfo(coloring,"Not calling MatZeroEntries()\n");CHKERRQ(ierr); 2129 } else { 2130 PetscTruth assembled; 2131 ierr = MatAssembled(J,&assembled);CHKERRQ(ierr); 2132 if (assembled) { 2133 ierr = MatZeroEntries(J);CHKERRQ(ierr); 2134 } 2135 } 2136 2137 ierr = VecGetOwnershipRange(x1,&start,&end);CHKERRQ(ierr); 2138 ierr = VecGetSize(x1,&N);CHKERRQ(ierr); 2139 2140 /* 2141 This is a horrible, horrible, hack. See DMMGComputeJacobian_Multigrid() it inproperly sets 2142 coloring->F for the coarser grids from the finest 2143 */ 2144 if (coloring->F) { 2145 ierr = VecGetLocalSize(coloring->F,&m1);CHKERRQ(ierr); 2146 ierr = VecGetLocalSize(w1,&m2);CHKERRQ(ierr); 2147 if (m1 != m2) { 2148 coloring->F = 0; 2149 } 2150 } 2151 2152 if (coloring->F) { 2153 w1 = coloring->F; 2154 coloring->F = 0; 2155 } else { 2156 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2157 ierr = (*f)(sctx,x1,w1,fctx);CHKERRQ(ierr); 2158 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2159 } 2160 2161 /* 2162 Compute all the scale factors and share with other processors 2163 */ 2164 ierr = VecGetArray(x1,&xx);CHKERRQ(ierr);xx = xx - start; 2165 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr);vscale_array = vscale_array - start; 2166 for (k=0; k<coloring->ncolors; k++) { 2167 /* 2168 Loop over each column associated with color adding the 2169 perturbation to the vector w3. 2170 */ 2171 for (l=0; l<coloring->ncolumns[k]; l++) { 2172 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2173 dx = xx[col]; 2174 if (dx == 0.0) dx = 1.0; 2175 #if !defined(PETSC_USE_COMPLEX) 2176 if (dx < umin && dx >= 0.0) dx = umin; 2177 else if (dx < 0.0 && dx > -umin) dx = -umin; 2178 #else 2179 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2180 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2181 #endif 2182 dx *= epsilon; 2183 vscale_array[col] = 1.0/dx; 2184 } 2185 } 2186 vscale_array = vscale_array + start;ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2187 ierr = VecGhostUpdateBegin(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2188 ierr = VecGhostUpdateEnd(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2189 2190 /* ierr = VecView(coloring->vscale,PETSC_VIEWER_STDOUT_WORLD); 2191 ierr = VecView(x1,PETSC_VIEWER_STDOUT_WORLD);*/ 2192 2193 if (coloring->vscaleforrow) vscaleforrow = coloring->vscaleforrow; 2194 else vscaleforrow = coloring->columnsforrow; 2195 2196 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2197 /* 2198 Loop over each color 2199 */ 2200 for (k=0; k<coloring->ncolors; k++) { 2201 coloring->currentcolor = k; 2202 ierr = VecCopy(x1,w3);CHKERRQ(ierr); 2203 ierr = VecGetArray(w3,&w3_array);CHKERRQ(ierr);w3_array = w3_array - start; 2204 /* 2205 Loop over each column associated with color adding the 2206 perturbation to the vector w3. 2207 */ 2208 for (l=0; l<coloring->ncolumns[k]; l++) { 2209 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2210 dx = xx[col]; 2211 if (dx == 0.0) dx = 1.0; 2212 #if !defined(PETSC_USE_COMPLEX) 2213 if (dx < umin && dx >= 0.0) dx = umin; 2214 else if (dx < 0.0 && dx > -umin) dx = -umin; 2215 #else 2216 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2217 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2218 #endif 2219 dx *= epsilon; 2220 if (!PetscAbsScalar(dx)) SETERRQ(PETSC_ERR_PLIB,"Computed 0 differencing parameter"); 2221 w3_array[col] += dx; 2222 } 2223 w3_array = w3_array + start; ierr = VecRestoreArray(w3,&w3_array);CHKERRQ(ierr); 2224 2225 /* 2226 Evaluate function at x1 + dx (here dx is a vector of perturbations) 2227 */ 2228 2229 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2230 ierr = (*f)(sctx,w3,w2,fctx);CHKERRQ(ierr); 2231 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2232 ierr = VecAXPY(w2,-1.0,w1);CHKERRQ(ierr); 2233 2234 /* 2235 Loop over rows of vector, putting results into Jacobian matrix 2236 */ 2237 ierr = VecGetArray(w2,&y);CHKERRQ(ierr); 2238 for (l=0; l<coloring->nrows[k]; l++) { 2239 row = coloring->rows[k][l]; 2240 col = coloring->columnsforrow[k][l]; 2241 y[row] *= vscale_array[vscaleforrow[k][l]]; 2242 srow = row + start; 2243 ierr = MatSetValues_SeqAIJ(J,1,&srow,1,&col,y+row,INSERT_VALUES);CHKERRQ(ierr); 2244 } 2245 ierr = VecRestoreArray(w2,&y);CHKERRQ(ierr); 2246 } 2247 coloring->currentcolor = k; 2248 ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2249 xx = xx + start; ierr = VecRestoreArray(x1,&xx);CHKERRQ(ierr); 2250 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2251 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2252 PetscFunctionReturn(0); 2253 } 2254 2255 #include "petscblaslapack.h" 2256 #undef __FUNCT__ 2257 #define __FUNCT__ "MatAXPY_SeqAIJ" 2258 PetscErrorCode MatAXPY_SeqAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 2259 { 2260 PetscErrorCode ierr; 2261 PetscInt i; 2262 Mat_SeqAIJ *x = (Mat_SeqAIJ *)X->data,*y = (Mat_SeqAIJ *)Y->data; 2263 PetscBLASInt one=1,bnz = PetscBLASIntCast(x->nz); 2264 2265 PetscFunctionBegin; 2266 if (str == SAME_NONZERO_PATTERN) { 2267 PetscScalar alpha = a; 2268 BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one); 2269 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 2270 if (y->xtoy && y->XtoY != X) { 2271 ierr = PetscFree(y->xtoy);CHKERRQ(ierr); 2272 ierr = MatDestroy(y->XtoY);CHKERRQ(ierr); 2273 } 2274 if (!y->xtoy) { /* get xtoy */ 2275 ierr = MatAXPYGetxtoy_Private(X->rmap->n,x->i,x->j,PETSC_NULL, y->i,y->j,PETSC_NULL, &y->xtoy);CHKERRQ(ierr); 2276 y->XtoY = X; 2277 ierr = PetscObjectReference((PetscObject)X);CHKERRQ(ierr); 2278 } 2279 for (i=0; i<x->nz; i++) y->a[y->xtoy[i]] += a*(x->a[i]); 2280 ierr = PetscInfo3(Y,"ratio of nnz(X)/nnz(Y): %d/%d = %G\n",x->nz,y->nz,(PetscReal)(x->nz)/y->nz);CHKERRQ(ierr); 2281 } else { 2282 ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr); 2283 } 2284 PetscFunctionReturn(0); 2285 } 2286 2287 #undef __FUNCT__ 2288 #define __FUNCT__ "MatSetBlockSize_SeqAIJ" 2289 PetscErrorCode MatSetBlockSize_SeqAIJ(Mat A,PetscInt bs) 2290 { 2291 PetscErrorCode ierr; 2292 2293 PetscFunctionBegin; 2294 ierr = PetscLayoutSetBlockSize(A->rmap,bs);CHKERRQ(ierr); 2295 ierr = PetscLayoutSetBlockSize(A->cmap,bs);CHKERRQ(ierr); 2296 PetscFunctionReturn(0); 2297 } 2298 2299 #undef __FUNCT__ 2300 #define __FUNCT__ "MatConjugate_SeqAIJ" 2301 PetscErrorCode PETSCMAT_DLLEXPORT MatConjugate_SeqAIJ(Mat mat) 2302 { 2303 #if defined(PETSC_USE_COMPLEX) 2304 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2305 PetscInt i,nz; 2306 PetscScalar *a; 2307 2308 PetscFunctionBegin; 2309 nz = aij->nz; 2310 a = aij->a; 2311 for (i=0; i<nz; i++) { 2312 a[i] = PetscConj(a[i]); 2313 } 2314 #else 2315 PetscFunctionBegin; 2316 #endif 2317 PetscFunctionReturn(0); 2318 } 2319 2320 #undef __FUNCT__ 2321 #define __FUNCT__ "MatGetRowMaxAbs_SeqAIJ" 2322 PetscErrorCode MatGetRowMaxAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2323 { 2324 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2325 PetscErrorCode ierr; 2326 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2327 PetscReal atmp; 2328 PetscScalar *x; 2329 MatScalar *aa; 2330 2331 PetscFunctionBegin; 2332 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2333 aa = a->a; 2334 ai = a->i; 2335 aj = a->j; 2336 2337 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2338 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2339 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2340 if (n != A->rmap->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2341 for (i=0; i<m; i++) { 2342 ncols = ai[1] - ai[0]; ai++; 2343 x[i] = 0.0; 2344 for (j=0; j<ncols; j++){ 2345 atmp = PetscAbsScalar(*aa); 2346 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2347 aa++; aj++; 2348 } 2349 } 2350 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2351 PetscFunctionReturn(0); 2352 } 2353 2354 #undef __FUNCT__ 2355 #define __FUNCT__ "MatGetRowMax_SeqAIJ" 2356 PetscErrorCode MatGetRowMax_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2357 { 2358 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2359 PetscErrorCode ierr; 2360 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2361 PetscScalar *x; 2362 MatScalar *aa; 2363 2364 PetscFunctionBegin; 2365 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2366 aa = a->a; 2367 ai = a->i; 2368 aj = a->j; 2369 2370 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2371 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2372 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2373 if (n != A->rmap->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2374 for (i=0; i<m; i++) { 2375 ncols = ai[1] - ai[0]; ai++; 2376 if (ncols == A->cmap->n) { /* row is dense */ 2377 x[i] = *aa; if (idx) idx[i] = 0; 2378 } else { /* row is sparse so already KNOW maximum is 0.0 or higher */ 2379 x[i] = 0.0; 2380 if (idx) { 2381 idx[i] = 0; /* in case ncols is zero */ 2382 for (j=0;j<ncols;j++) { /* find first implicit 0.0 in the row */ 2383 if (aj[j] > j) { 2384 idx[i] = j; 2385 break; 2386 } 2387 } 2388 } 2389 } 2390 for (j=0; j<ncols; j++){ 2391 if (PetscRealPart(x[i]) < PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2392 aa++; aj++; 2393 } 2394 } 2395 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2396 PetscFunctionReturn(0); 2397 } 2398 2399 #undef __FUNCT__ 2400 #define __FUNCT__ "MatGetRowMinAbs_SeqAIJ" 2401 PetscErrorCode MatGetRowMinAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2402 { 2403 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2404 PetscErrorCode ierr; 2405 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2406 PetscReal atmp; 2407 PetscScalar *x; 2408 MatScalar *aa; 2409 2410 PetscFunctionBegin; 2411 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2412 aa = a->a; 2413 ai = a->i; 2414 aj = a->j; 2415 2416 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2417 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2418 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2419 if (n != A->rmap->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2420 for (i=0; i<m; i++) { 2421 ncols = ai[1] - ai[0]; ai++; 2422 if (ncols) { 2423 /* Get first nonzero */ 2424 for(j = 0; j < ncols; j++) { 2425 atmp = PetscAbsScalar(aa[j]); 2426 if (atmp > 1.0e-12) {x[i] = atmp; if (idx) idx[i] = aj[j]; break;} 2427 } 2428 if (j == ncols) {x[i] = *aa; if (idx) idx[i] = *aj;} 2429 } else { 2430 x[i] = 0.0; if (idx) idx[i] = 0; 2431 } 2432 for(j = 0; j < ncols; j++) { 2433 atmp = PetscAbsScalar(*aa); 2434 if (atmp > 1.0e-12 && PetscAbsScalar(x[i]) > atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2435 aa++; aj++; 2436 } 2437 } 2438 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2439 PetscFunctionReturn(0); 2440 } 2441 2442 #undef __FUNCT__ 2443 #define __FUNCT__ "MatGetRowMin_SeqAIJ" 2444 PetscErrorCode MatGetRowMin_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2445 { 2446 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2447 PetscErrorCode ierr; 2448 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2449 PetscScalar *x; 2450 MatScalar *aa; 2451 2452 PetscFunctionBegin; 2453 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2454 aa = a->a; 2455 ai = a->i; 2456 aj = a->j; 2457 2458 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2459 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2460 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2461 if (n != A->rmap->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2462 for (i=0; i<m; i++) { 2463 ncols = ai[1] - ai[0]; ai++; 2464 if (ncols == A->cmap->n) { /* row is dense */ 2465 x[i] = *aa; if (idx) idx[i] = 0; 2466 } else { /* row is sparse so already KNOW minimum is 0.0 or lower */ 2467 x[i] = 0.0; 2468 if (idx) { /* find first implicit 0.0 in the row */ 2469 idx[i] = 0; /* in case ncols is zero */ 2470 for (j=0;j<ncols;j++) { 2471 if (aj[j] > j) { 2472 idx[i] = j; 2473 break; 2474 } 2475 } 2476 } 2477 } 2478 for (j=0; j<ncols; j++){ 2479 if (PetscRealPart(x[i]) > PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2480 aa++; aj++; 2481 } 2482 } 2483 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2484 PetscFunctionReturn(0); 2485 } 2486 extern PetscErrorCode PETSCMAT_DLLEXPORT MatFDColoringApply_AIJ(Mat,MatFDColoring,Vec,MatStructure*,void*); 2487 /* -------------------------------------------------------------------*/ 2488 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 2489 MatGetRow_SeqAIJ, 2490 MatRestoreRow_SeqAIJ, 2491 MatMult_SeqAIJ, 2492 /* 4*/ MatMultAdd_SeqAIJ, 2493 MatMultTranspose_SeqAIJ, 2494 MatMultTransposeAdd_SeqAIJ, 2495 0, 2496 0, 2497 0, 2498 /*10*/ 0, 2499 MatLUFactor_SeqAIJ, 2500 0, 2501 MatSOR_SeqAIJ, 2502 MatTranspose_SeqAIJ, 2503 /*15*/ MatGetInfo_SeqAIJ, 2504 MatEqual_SeqAIJ, 2505 MatGetDiagonal_SeqAIJ, 2506 MatDiagonalScale_SeqAIJ, 2507 MatNorm_SeqAIJ, 2508 /*20*/ 0, 2509 MatAssemblyEnd_SeqAIJ, 2510 MatSetOption_SeqAIJ, 2511 MatZeroEntries_SeqAIJ, 2512 /*24*/ MatZeroRows_SeqAIJ, 2513 0, 2514 0, 2515 0, 2516 0, 2517 /*29*/ MatSetUpPreallocation_SeqAIJ, 2518 0, 2519 0, 2520 MatGetArray_SeqAIJ, 2521 MatRestoreArray_SeqAIJ, 2522 /*34*/ MatDuplicate_SeqAIJ, 2523 0, 2524 0, 2525 MatILUFactor_SeqAIJ, 2526 0, 2527 /*39*/ MatAXPY_SeqAIJ, 2528 MatGetSubMatrices_SeqAIJ, 2529 MatIncreaseOverlap_SeqAIJ, 2530 MatGetValues_SeqAIJ, 2531 MatCopy_SeqAIJ, 2532 /*44*/ MatGetRowMax_SeqAIJ, 2533 MatScale_SeqAIJ, 2534 0, 2535 MatDiagonalSet_SeqAIJ, 2536 MatILUDTFactor_SeqAIJ, 2537 /*49*/ MatSetBlockSize_SeqAIJ, 2538 MatGetRowIJ_SeqAIJ, 2539 MatRestoreRowIJ_SeqAIJ, 2540 MatGetColumnIJ_SeqAIJ, 2541 MatRestoreColumnIJ_SeqAIJ, 2542 /*54*/ MatFDColoringCreate_SeqAIJ, 2543 0, 2544 0, 2545 MatPermute_SeqAIJ, 2546 0, 2547 /*59*/ 0, 2548 MatDestroy_SeqAIJ, 2549 MatView_SeqAIJ, 2550 0, 2551 0, 2552 /*64*/ 0, 2553 0, 2554 0, 2555 0, 2556 0, 2557 /*69*/ MatGetRowMaxAbs_SeqAIJ, 2558 MatGetRowMinAbs_SeqAIJ, 2559 0, 2560 MatSetColoring_SeqAIJ, 2561 #if defined(PETSC_HAVE_ADIC) 2562 MatSetValuesAdic_SeqAIJ, 2563 #else 2564 0, 2565 #endif 2566 /*74*/ MatSetValuesAdifor_SeqAIJ, 2567 MatFDColoringApply_AIJ, 2568 0, 2569 0, 2570 0, 2571 /*79*/ 0, 2572 0, 2573 0, 2574 0, 2575 MatLoad_SeqAIJ, 2576 /*84*/ MatIsSymmetric_SeqAIJ, 2577 MatIsHermitian_SeqAIJ, 2578 0, 2579 0, 2580 0, 2581 /*89*/ MatMatMult_SeqAIJ_SeqAIJ, 2582 MatMatMultSymbolic_SeqAIJ_SeqAIJ, 2583 MatMatMultNumeric_SeqAIJ_SeqAIJ, 2584 MatPtAP_Basic, 2585 MatPtAPSymbolic_SeqAIJ, 2586 /*94*/ MatPtAPNumeric_SeqAIJ, 2587 MatMatMultTranspose_SeqAIJ_SeqAIJ, 2588 MatMatMultTransposeSymbolic_SeqAIJ_SeqAIJ, 2589 MatMatMultTransposeNumeric_SeqAIJ_SeqAIJ, 2590 MatPtAPSymbolic_SeqAIJ_SeqAIJ, 2591 /*99*/ MatPtAPNumeric_SeqAIJ_SeqAIJ, 2592 0, 2593 0, 2594 MatConjugate_SeqAIJ, 2595 0, 2596 /*104*/MatSetValuesRow_SeqAIJ, 2597 MatRealPart_SeqAIJ, 2598 MatImaginaryPart_SeqAIJ, 2599 0, 2600 0, 2601 /*109*/0, 2602 0, 2603 MatGetRowMin_SeqAIJ, 2604 0, 2605 MatMissingDiagonal_SeqAIJ, 2606 /*114*/0, 2607 0, 2608 0, 2609 MatILUDTFactorSymbolic_SeqAIJ, 2610 MatILUDTFactorNumeric_SeqAIJ 2611 }; 2612 2613 EXTERN_C_BEGIN 2614 #undef __FUNCT__ 2615 #define __FUNCT__ "MatSeqAIJSetColumnIndices_SeqAIJ" 2616 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 2617 { 2618 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2619 PetscInt i,nz,n; 2620 2621 PetscFunctionBegin; 2622 2623 nz = aij->maxnz; 2624 n = mat->rmap->n; 2625 for (i=0; i<nz; i++) { 2626 aij->j[i] = indices[i]; 2627 } 2628 aij->nz = nz; 2629 for (i=0; i<n; i++) { 2630 aij->ilen[i] = aij->imax[i]; 2631 } 2632 2633 PetscFunctionReturn(0); 2634 } 2635 EXTERN_C_END 2636 2637 #undef __FUNCT__ 2638 #define __FUNCT__ "MatSeqAIJSetColumnIndices" 2639 /*@ 2640 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 2641 in the matrix. 2642 2643 Input Parameters: 2644 + mat - the SeqAIJ matrix 2645 - indices - the column indices 2646 2647 Level: advanced 2648 2649 Notes: 2650 This can be called if you have precomputed the nonzero structure of the 2651 matrix and want to provide it to the matrix object to improve the performance 2652 of the MatSetValues() operation. 2653 2654 You MUST have set the correct numbers of nonzeros per row in the call to 2655 MatCreateSeqAIJ(), and the columns indices MUST be sorted. 2656 2657 MUST be called before any calls to MatSetValues(); 2658 2659 The indices should start with zero, not one. 2660 2661 @*/ 2662 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 2663 { 2664 PetscErrorCode ierr,(*f)(Mat,PetscInt *); 2665 2666 PetscFunctionBegin; 2667 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 2668 PetscValidPointer(indices,2); 2669 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqAIJSetColumnIndices_C",(void (**)(void))&f);CHKERRQ(ierr); 2670 if (f) { 2671 ierr = (*f)(mat,indices);CHKERRQ(ierr); 2672 } else { 2673 SETERRQ(PETSC_ERR_SUP,"Wrong type of matrix to set column indices"); 2674 } 2675 PetscFunctionReturn(0); 2676 } 2677 2678 /* ----------------------------------------------------------------------------------------*/ 2679 2680 EXTERN_C_BEGIN 2681 #undef __FUNCT__ 2682 #define __FUNCT__ "MatStoreValues_SeqAIJ" 2683 PetscErrorCode PETSCMAT_DLLEXPORT MatStoreValues_SeqAIJ(Mat mat) 2684 { 2685 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2686 PetscErrorCode ierr; 2687 size_t nz = aij->i[mat->rmap->n]; 2688 2689 PetscFunctionBegin; 2690 if (aij->nonew != 1) { 2691 SETERRQ(PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 2692 } 2693 2694 /* allocate space for values if not already there */ 2695 if (!aij->saved_values) { 2696 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 2697 ierr = PetscLogObjectMemory(mat,(nz+1)*sizeof(PetscScalar));CHKERRQ(ierr); 2698 } 2699 2700 /* copy values over */ 2701 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2702 PetscFunctionReturn(0); 2703 } 2704 EXTERN_C_END 2705 2706 #undef __FUNCT__ 2707 #define __FUNCT__ "MatStoreValues" 2708 /*@ 2709 MatStoreValues - Stashes a copy of the matrix values; this allows, for 2710 example, reuse of the linear part of a Jacobian, while recomputing the 2711 nonlinear portion. 2712 2713 Collect on Mat 2714 2715 Input Parameters: 2716 . mat - the matrix (currently only AIJ matrices support this option) 2717 2718 Level: advanced 2719 2720 Common Usage, with SNESSolve(): 2721 $ Create Jacobian matrix 2722 $ Set linear terms into matrix 2723 $ Apply boundary conditions to matrix, at this time matrix must have 2724 $ final nonzero structure (i.e. setting the nonlinear terms and applying 2725 $ boundary conditions again will not change the nonzero structure 2726 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 2727 $ ierr = MatStoreValues(mat); 2728 $ Call SNESSetJacobian() with matrix 2729 $ In your Jacobian routine 2730 $ ierr = MatRetrieveValues(mat); 2731 $ Set nonlinear terms in matrix 2732 2733 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 2734 $ // build linear portion of Jacobian 2735 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 2736 $ ierr = MatStoreValues(mat); 2737 $ loop over nonlinear iterations 2738 $ ierr = MatRetrieveValues(mat); 2739 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 2740 $ // call MatAssemblyBegin/End() on matrix 2741 $ Solve linear system with Jacobian 2742 $ endloop 2743 2744 Notes: 2745 Matrix must already be assemblied before calling this routine 2746 Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before 2747 calling this routine. 2748 2749 When this is called multiple times it overwrites the previous set of stored values 2750 and does not allocated additional space. 2751 2752 .seealso: MatRetrieveValues() 2753 2754 @*/ 2755 PetscErrorCode PETSCMAT_DLLEXPORT MatStoreValues(Mat mat) 2756 { 2757 PetscErrorCode ierr,(*f)(Mat); 2758 2759 PetscFunctionBegin; 2760 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 2761 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2762 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2763 2764 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatStoreValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2765 if (f) { 2766 ierr = (*f)(mat);CHKERRQ(ierr); 2767 } else { 2768 SETERRQ(PETSC_ERR_SUP,"Wrong type of matrix to store values"); 2769 } 2770 PetscFunctionReturn(0); 2771 } 2772 2773 EXTERN_C_BEGIN 2774 #undef __FUNCT__ 2775 #define __FUNCT__ "MatRetrieveValues_SeqAIJ" 2776 PetscErrorCode PETSCMAT_DLLEXPORT MatRetrieveValues_SeqAIJ(Mat mat) 2777 { 2778 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2779 PetscErrorCode ierr; 2780 PetscInt nz = aij->i[mat->rmap->n]; 2781 2782 PetscFunctionBegin; 2783 if (aij->nonew != 1) { 2784 SETERRQ(PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 2785 } 2786 if (!aij->saved_values) { 2787 SETERRQ(PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 2788 } 2789 /* copy values over */ 2790 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2791 PetscFunctionReturn(0); 2792 } 2793 EXTERN_C_END 2794 2795 #undef __FUNCT__ 2796 #define __FUNCT__ "MatRetrieveValues" 2797 /*@ 2798 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 2799 example, reuse of the linear part of a Jacobian, while recomputing the 2800 nonlinear portion. 2801 2802 Collect on Mat 2803 2804 Input Parameters: 2805 . mat - the matrix (currently on AIJ matrices support this option) 2806 2807 Level: advanced 2808 2809 .seealso: MatStoreValues() 2810 2811 @*/ 2812 PetscErrorCode PETSCMAT_DLLEXPORT MatRetrieveValues(Mat mat) 2813 { 2814 PetscErrorCode ierr,(*f)(Mat); 2815 2816 PetscFunctionBegin; 2817 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 2818 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2819 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2820 2821 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatRetrieveValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2822 if (f) { 2823 ierr = (*f)(mat);CHKERRQ(ierr); 2824 } else { 2825 SETERRQ(PETSC_ERR_SUP,"Wrong type of matrix to retrieve values"); 2826 } 2827 PetscFunctionReturn(0); 2828 } 2829 2830 2831 /* --------------------------------------------------------------------------------*/ 2832 #undef __FUNCT__ 2833 #define __FUNCT__ "MatCreateSeqAIJ" 2834 /*@C 2835 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 2836 (the default parallel PETSc format). For good matrix assembly performance 2837 the user should preallocate the matrix storage by setting the parameter nz 2838 (or the array nnz). By setting these parameters accurately, performance 2839 during matrix assembly can be increased by more than a factor of 50. 2840 2841 Collective on MPI_Comm 2842 2843 Input Parameters: 2844 + comm - MPI communicator, set to PETSC_COMM_SELF 2845 . m - number of rows 2846 . n - number of columns 2847 . nz - number of nonzeros per row (same for all rows) 2848 - nnz - array containing the number of nonzeros in the various rows 2849 (possibly different for each row) or PETSC_NULL 2850 2851 Output Parameter: 2852 . A - the matrix 2853 2854 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 2855 MatXXXXSetPreallocation() paradgm instead of this routine directly. 2856 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 2857 2858 Notes: 2859 If nnz is given then nz is ignored 2860 2861 The AIJ format (also called the Yale sparse matrix format or 2862 compressed row storage), is fully compatible with standard Fortran 77 2863 storage. That is, the stored row and column indices can begin at 2864 either one (as in Fortran) or zero. See the users' manual for details. 2865 2866 Specify the preallocated storage with either nz or nnz (not both). 2867 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2868 allocation. For large problems you MUST preallocate memory or you 2869 will get TERRIBLE performance, see the users' manual chapter on matrices. 2870 2871 By default, this format uses inodes (identical nodes) when possible, to 2872 improve numerical efficiency of matrix-vector products and solves. We 2873 search for consecutive rows with the same nonzero structure, thereby 2874 reusing matrix information to achieve increased efficiency. 2875 2876 Options Database Keys: 2877 + -mat_no_inode - Do not use inodes 2878 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 2879 - -mat_aij_oneindex - Internally use indexing starting at 1 2880 rather than 0. Note that when calling MatSetValues(), 2881 the user still MUST index entries starting at 0! 2882 2883 Level: intermediate 2884 2885 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2886 2887 @*/ 2888 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 2889 { 2890 PetscErrorCode ierr; 2891 2892 PetscFunctionBegin; 2893 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2894 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2895 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 2896 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,(PetscInt*)nnz);CHKERRQ(ierr); 2897 PetscFunctionReturn(0); 2898 } 2899 2900 #undef __FUNCT__ 2901 #define __FUNCT__ "MatSeqAIJSetPreallocation" 2902 /*@C 2903 MatSeqAIJSetPreallocation - For good matrix assembly performance 2904 the user should preallocate the matrix storage by setting the parameter nz 2905 (or the array nnz). By setting these parameters accurately, performance 2906 during matrix assembly can be increased by more than a factor of 50. 2907 2908 Collective on MPI_Comm 2909 2910 Input Parameters: 2911 + B - The matrix-free 2912 . nz - number of nonzeros per row (same for all rows) 2913 - nnz - array containing the number of nonzeros in the various rows 2914 (possibly different for each row) or PETSC_NULL 2915 2916 Notes: 2917 If nnz is given then nz is ignored 2918 2919 The AIJ format (also called the Yale sparse matrix format or 2920 compressed row storage), is fully compatible with standard Fortran 77 2921 storage. That is, the stored row and column indices can begin at 2922 either one (as in Fortran) or zero. See the users' manual for details. 2923 2924 Specify the preallocated storage with either nz or nnz (not both). 2925 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2926 allocation. For large problems you MUST preallocate memory or you 2927 will get TERRIBLE performance, see the users' manual chapter on matrices. 2928 2929 You can call MatGetInfo() to get information on how effective the preallocation was; 2930 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 2931 You can also run with the option -info and look for messages with the string 2932 malloc in them to see if additional memory allocation was needed. 2933 2934 Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix 2935 entries or columns indices 2936 2937 By default, this format uses inodes (identical nodes) when possible, to 2938 improve numerical efficiency of matrix-vector products and solves. We 2939 search for consecutive rows with the same nonzero structure, thereby 2940 reusing matrix information to achieve increased efficiency. 2941 2942 Options Database Keys: 2943 + -mat_no_inode - Do not use inodes 2944 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 2945 - -mat_aij_oneindex - Internally use indexing starting at 1 2946 rather than 0. Note that when calling MatSetValues(), 2947 the user still MUST index entries starting at 0! 2948 2949 Level: intermediate 2950 2951 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatGetInfo() 2952 2953 @*/ 2954 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 2955 { 2956 PetscErrorCode ierr,(*f)(Mat,PetscInt,const PetscInt[]); 2957 2958 PetscFunctionBegin; 2959 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 2960 if (f) { 2961 ierr = (*f)(B,nz,nnz);CHKERRQ(ierr); 2962 } 2963 PetscFunctionReturn(0); 2964 } 2965 2966 EXTERN_C_BEGIN 2967 #undef __FUNCT__ 2968 #define __FUNCT__ "MatSeqAIJSetPreallocation_SeqAIJ" 2969 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,PetscInt *nnz) 2970 { 2971 Mat_SeqAIJ *b; 2972 PetscTruth skipallocation = PETSC_FALSE; 2973 PetscErrorCode ierr; 2974 PetscInt i; 2975 2976 PetscFunctionBegin; 2977 2978 if (nz == MAT_SKIP_ALLOCATION) { 2979 skipallocation = PETSC_TRUE; 2980 nz = 0; 2981 } 2982 2983 ierr = PetscLayoutSetBlockSize(B->rmap,1);CHKERRQ(ierr); 2984 ierr = PetscLayoutSetBlockSize(B->cmap,1);CHKERRQ(ierr); 2985 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 2986 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 2987 2988 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 2989 if (nz < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 2990 if (nnz) { 2991 for (i=0; i<B->rmap->n; i++) { 2992 if (nnz[i] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]); 2993 if (nnz[i] > B->cmap->n) SETERRQ3(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be greater than row length: local row %d value %d rowlength %d",i,nnz[i],B->cmap->n); 2994 } 2995 } 2996 2997 B->preallocated = PETSC_TRUE; 2998 b = (Mat_SeqAIJ*)B->data; 2999 3000 if (!skipallocation) { 3001 if (!b->imax) { 3002 ierr = PetscMalloc2(B->rmap->n,PetscInt,&b->imax,B->rmap->n,PetscInt,&b->ilen);CHKERRQ(ierr); 3003 ierr = PetscLogObjectMemory(B,2*B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3004 } 3005 if (!nnz) { 3006 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 3007 else if (nz <= 0) nz = 1; 3008 for (i=0; i<B->rmap->n; i++) b->imax[i] = nz; 3009 nz = nz*B->rmap->n; 3010 } else { 3011 nz = 0; 3012 for (i=0; i<B->rmap->n; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 3013 } 3014 /* b->ilen will count nonzeros in each row so far. */ 3015 for (i=0; i<B->rmap->n; i++) { b->ilen[i] = 0; } 3016 3017 /* allocate the matrix space */ 3018 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 3019 ierr = PetscMalloc3(nz,PetscScalar,&b->a,nz,PetscInt,&b->j,B->rmap->n+1,PetscInt,&b->i);CHKERRQ(ierr); 3020 ierr = PetscLogObjectMemory(B,(B->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 3021 b->i[0] = 0; 3022 for (i=1; i<B->rmap->n+1; i++) { 3023 b->i[i] = b->i[i-1] + b->imax[i-1]; 3024 } 3025 b->singlemalloc = PETSC_TRUE; 3026 b->free_a = PETSC_TRUE; 3027 b->free_ij = PETSC_TRUE; 3028 } else { 3029 b->free_a = PETSC_FALSE; 3030 b->free_ij = PETSC_FALSE; 3031 } 3032 3033 b->nz = 0; 3034 b->maxnz = nz; 3035 B->info.nz_unneeded = (double)b->maxnz; 3036 PetscFunctionReturn(0); 3037 } 3038 EXTERN_C_END 3039 3040 #undef __FUNCT__ 3041 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR" 3042 /*@ 3043 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 3044 3045 Input Parameters: 3046 + B - the matrix 3047 . i - the indices into j for the start of each row (starts with zero) 3048 . j - the column indices for each row (starts with zero) these must be sorted for each row 3049 - v - optional values in the matrix 3050 3051 Level: developer 3052 3053 The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays() 3054 3055 .keywords: matrix, aij, compressed row, sparse, sequential 3056 3057 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), SeqAIJ 3058 @*/ 3059 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 3060 { 3061 PetscErrorCode (*f)(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]); 3062 PetscErrorCode ierr; 3063 3064 PetscFunctionBegin; 3065 PetscValidHeaderSpecific(B,MAT_COOKIE,1); 3066 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C",(void (**)(void))&f);CHKERRQ(ierr); 3067 if (f) { 3068 ierr = (*f)(B,i,j,v);CHKERRQ(ierr); 3069 } 3070 PetscFunctionReturn(0); 3071 } 3072 3073 EXTERN_C_BEGIN 3074 #undef __FUNCT__ 3075 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR_SeqAIJ" 3076 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[]) 3077 { 3078 PetscInt i; 3079 PetscInt m,n; 3080 PetscInt nz; 3081 PetscInt *nnz, nz_max = 0; 3082 PetscScalar *values; 3083 PetscErrorCode ierr; 3084 3085 PetscFunctionBegin; 3086 ierr = MatGetSize(B, &m, &n);CHKERRQ(ierr); 3087 3088 if (Ii[0]) { 3089 SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %D", Ii[0]); 3090 } 3091 ierr = PetscMalloc((m+1) * sizeof(PetscInt), &nnz);CHKERRQ(ierr); 3092 for(i = 0; i < m; i++) { 3093 nz = Ii[i+1]- Ii[i]; 3094 nz_max = PetscMax(nz_max, nz); 3095 if (nz < 0) { 3096 SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Local row %D has a negative number of columns %D", i, nnz); 3097 } 3098 nnz[i] = nz; 3099 } 3100 ierr = MatSeqAIJSetPreallocation(B, 0, nnz);CHKERRQ(ierr); 3101 ierr = PetscFree(nnz);CHKERRQ(ierr); 3102 3103 if (v) { 3104 values = (PetscScalar*) v; 3105 } else { 3106 ierr = PetscMalloc(nz_max*sizeof(PetscScalar), &values);CHKERRQ(ierr); 3107 ierr = PetscMemzero(values, nz_max*sizeof(PetscScalar));CHKERRQ(ierr); 3108 } 3109 3110 for(i = 0; i < m; i++) { 3111 nz = Ii[i+1] - Ii[i]; 3112 ierr = MatSetValues_SeqAIJ(B, 1, &i, nz, J+Ii[i], values + (v ? Ii[i] : 0), INSERT_VALUES);CHKERRQ(ierr); 3113 } 3114 3115 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3116 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3117 3118 if (!v) { 3119 ierr = PetscFree(values);CHKERRQ(ierr); 3120 } 3121 PetscFunctionReturn(0); 3122 } 3123 EXTERN_C_END 3124 3125 #include "../src/mat/impls/dense/seq/dense.h" 3126 #include "private/petscaxpy.h" 3127 3128 #undef __FUNCT__ 3129 #define __FUNCT__ "MatMatMultNumeric_SeqDense_SeqAIJ" 3130 /* 3131 Computes (B'*A')' since computing B*A directly is untenable 3132 3133 n p p 3134 ( ) ( ) ( ) 3135 m ( A ) * n ( B ) = m ( C ) 3136 ( ) ( ) ( ) 3137 3138 */ 3139 PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C) 3140 { 3141 PetscErrorCode ierr; 3142 Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data; 3143 Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data; 3144 Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data; 3145 PetscInt i,n,m,q,p; 3146 const PetscInt *ii,*idx; 3147 const PetscScalar *b,*a,*a_q; 3148 PetscScalar *c,*c_q; 3149 3150 PetscFunctionBegin; 3151 m = A->rmap->n; 3152 n = A->cmap->n; 3153 p = B->cmap->n; 3154 a = sub_a->v; 3155 b = sub_b->a; 3156 c = sub_c->v; 3157 ierr = PetscMemzero(c,m*p*sizeof(PetscScalar));CHKERRQ(ierr); 3158 3159 ii = sub_b->i; 3160 idx = sub_b->j; 3161 for (i=0; i<n; i++) { 3162 q = ii[i+1] - ii[i]; 3163 while (q-->0) { 3164 c_q = c + m*(*idx); 3165 a_q = a + m*i; 3166 PetscAXPY(c_q,*b,a_q,m); 3167 idx++; 3168 b++; 3169 } 3170 } 3171 PetscFunctionReturn(0); 3172 } 3173 3174 #undef __FUNCT__ 3175 #define __FUNCT__ "MatMatMultSymbolic_SeqDense_SeqAIJ" 3176 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 3177 { 3178 PetscErrorCode ierr; 3179 PetscInt m=A->rmap->n,n=B->cmap->n; 3180 Mat Cmat; 3181 3182 PetscFunctionBegin; 3183 if (A->cmap->n != B->rmap->n) SETERRQ2(PETSC_ERR_ARG_SIZ,"A->cmap->n %d != B->rmap->n %d\n",A->cmap->n,B->rmap->n); 3184 ierr = MatCreate(((PetscObject)A)->comm,&Cmat);CHKERRQ(ierr); 3185 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 3186 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 3187 ierr = MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);CHKERRQ(ierr); 3188 Cmat->assembled = PETSC_TRUE; 3189 *C = Cmat; 3190 PetscFunctionReturn(0); 3191 } 3192 3193 /* ----------------------------------------------------------------*/ 3194 #undef __FUNCT__ 3195 #define __FUNCT__ "MatMatMult_SeqDense_SeqAIJ" 3196 PetscErrorCode MatMatMult_SeqDense_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 3197 { 3198 PetscErrorCode ierr; 3199 3200 PetscFunctionBegin; 3201 if (scall == MAT_INITIAL_MATRIX){ 3202 ierr = MatMatMultSymbolic_SeqDense_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 3203 } 3204 ierr = MatMatMultNumeric_SeqDense_SeqAIJ(A,B,*C);CHKERRQ(ierr); 3205 PetscFunctionReturn(0); 3206 } 3207 3208 3209 /*MC 3210 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 3211 based on compressed sparse row format. 3212 3213 Options Database Keys: 3214 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 3215 3216 Level: beginner 3217 3218 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 3219 M*/ 3220 3221 EXTERN_C_BEGIN 3222 #if defined(PETSC_HAVE_PASTIX) 3223 extern PetscErrorCode MatGetFactor_seqaij_pastix(Mat,MatFactorType,Mat*); 3224 #endif 3225 #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SCALAR_SINGLE) && !defined(PETSC_USE_SCALAR_MAT_SINGLE) 3226 extern PetscErrorCode MatGetFactor_seqaij_essl(Mat,MatFactorType,Mat *); 3227 #endif 3228 extern PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqAIJ_SeqCRL(Mat,MatType,MatReuse,Mat*); 3229 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_petsc(Mat,MatFactorType,Mat*); 3230 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_bas(Mat,MatFactorType,Mat*); 3231 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactorAvailable_seqaij_petsc(Mat,MatFactorType,PetscTruth *); 3232 #if defined(PETSC_HAVE_MUMPS) 3233 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_mumps(Mat,MatFactorType,Mat*); 3234 #endif 3235 #if defined(PETSC_HAVE_SUPERLU) 3236 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_superlu(Mat,MatFactorType,Mat*); 3237 #endif 3238 #if defined(PETSC_HAVE_SUPERLU_DIST) 3239 extern PetscErrorCode MatGetFactor_seqaij_superlu_dist(Mat,MatFactorType,Mat*); 3240 #endif 3241 #if defined(PETSC_HAVE_SPOOLES) 3242 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_spooles(Mat,MatFactorType,Mat*); 3243 #endif 3244 #if defined(PETSC_HAVE_UMFPACK) 3245 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_umfpack(Mat,MatFactorType,Mat*); 3246 #endif 3247 #if defined(PETSC_HAVE_LUSOL) 3248 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_lusol(Mat,MatFactorType,Mat*); 3249 #endif 3250 EXTERN_C_END 3251 3252 3253 EXTERN_C_BEGIN 3254 #undef __FUNCT__ 3255 #define __FUNCT__ "MatCreate_SeqAIJ" 3256 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_SeqAIJ(Mat B) 3257 { 3258 Mat_SeqAIJ *b; 3259 PetscErrorCode ierr; 3260 PetscMPIInt size; 3261 3262 PetscFunctionBegin; 3263 ierr = MPI_Comm_size(((PetscObject)B)->comm,&size);CHKERRQ(ierr); 3264 if (size > 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 3265 3266 ierr = PetscNewLog(B,Mat_SeqAIJ,&b);CHKERRQ(ierr); 3267 B->data = (void*)b; 3268 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 3269 B->mapping = 0; 3270 b->row = 0; 3271 b->col = 0; 3272 b->icol = 0; 3273 b->reallocs = 0; 3274 b->ignorezeroentries = PETSC_FALSE; 3275 b->roworiented = PETSC_TRUE; 3276 b->nonew = 0; 3277 b->diag = 0; 3278 b->solve_work = 0; 3279 B->spptr = 0; 3280 b->saved_values = 0; 3281 b->idiag = 0; 3282 b->mdiag = 0; 3283 b->ssor_work = 0; 3284 b->omega = 1.0; 3285 b->fshift = 0.0; 3286 b->idiagvalid = PETSC_FALSE; 3287 b->keepnonzeropattern = PETSC_FALSE; 3288 b->xtoy = 0; 3289 b->XtoY = 0; 3290 b->compressedrow.use = PETSC_FALSE; 3291 b->compressedrow.nrows = B->rmap->n; 3292 b->compressedrow.i = PETSC_NULL; 3293 b->compressedrow.rindex = PETSC_NULL; 3294 b->compressedrow.checked = PETSC_FALSE; 3295 B->same_nonzero = PETSC_FALSE; 3296 3297 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3298 #if defined(PETSC_HAVE_PASTIX) 3299 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_pastix_C", 3300 "MatGetFactor_seqaij_pastix", 3301 MatGetFactor_seqaij_pastix);CHKERRQ(ierr); 3302 #endif 3303 #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SCALAR_SINGLE) && !defined(PETSC_USE_SCALAR_MAT_SINGLE) 3304 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_essl_C", 3305 "MatGetFactor_seqaij_essl", 3306 MatGetFactor_seqaij_essl);CHKERRQ(ierr); 3307 #endif 3308 #if defined(PETSC_HAVE_SUPERLU) 3309 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_superlu_C", 3310 "MatGetFactor_seqaij_superlu", 3311 MatGetFactor_seqaij_superlu);CHKERRQ(ierr); 3312 #endif 3313 #if defined(PETSC_HAVE_SUPERLU_DIST) 3314 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_superlu_dist_C", 3315 "MatGetFactor_seqaij_superlu_dist", 3316 MatGetFactor_seqaij_superlu_dist);CHKERRQ(ierr); 3317 #endif 3318 #if defined(PETSC_HAVE_SPOOLES) 3319 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_spooles_C", 3320 "MatGetFactor_seqaij_spooles", 3321 MatGetFactor_seqaij_spooles);CHKERRQ(ierr); 3322 #endif 3323 #if defined(PETSC_HAVE_MUMPS) 3324 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_mumps_C", 3325 "MatGetFactor_seqaij_mumps", 3326 MatGetFactor_seqaij_mumps);CHKERRQ(ierr); 3327 #endif 3328 #if defined(PETSC_HAVE_UMFPACK) 3329 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_umfpack_C", 3330 "MatGetFactor_seqaij_umfpack", 3331 MatGetFactor_seqaij_umfpack);CHKERRQ(ierr); 3332 #endif 3333 #if defined(PETSC_HAVE_LUSOL) 3334 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_lusol_C", 3335 "MatGetFactor_seqaij_lusol", 3336 MatGetFactor_seqaij_lusol);CHKERRQ(ierr); 3337 #endif 3338 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_petsc_C", 3339 "MatGetFactor_seqaij_petsc", 3340 MatGetFactor_seqaij_petsc);CHKERRQ(ierr); 3341 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactorAvailable_petsc_C", 3342 "MatGetFactorAvailable_seqaij_petsc", 3343 MatGetFactorAvailable_seqaij_petsc);CHKERRQ(ierr); 3344 /* ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_bas_C", 3345 "MatGetFactor_seqaij_bas", 3346 MatGetFactor_seqaij_bas); */ 3347 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 3348 "MatSeqAIJSetColumnIndices_SeqAIJ", 3349 MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 3350 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 3351 "MatStoreValues_SeqAIJ", 3352 MatStoreValues_SeqAIJ);CHKERRQ(ierr); 3353 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 3354 "MatRetrieveValues_SeqAIJ", 3355 MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 3356 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqsbaij_C", 3357 "MatConvert_SeqAIJ_SeqSBAIJ", 3358 MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 3359 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqbaij_C", 3360 "MatConvert_SeqAIJ_SeqBAIJ", 3361 MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 3362 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqcsrperm_C", 3363 "MatConvert_SeqAIJ_SeqCSRPERM", 3364 MatConvert_SeqAIJ_SeqCSRPERM);CHKERRQ(ierr); 3365 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqcrl_C", 3366 "MatConvert_SeqAIJ_SeqCRL", 3367 MatConvert_SeqAIJ_SeqCRL);CHKERRQ(ierr); 3368 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsTranspose_C", 3369 "MatIsTranspose_SeqAIJ", 3370 MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3371 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsHermitianTranspose_C", 3372 "MatIsHermitianTranspose_SeqAIJ", 3373 MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3374 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocation_C", 3375 "MatSeqAIJSetPreallocation_SeqAIJ", 3376 MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 3377 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C", 3378 "MatSeqAIJSetPreallocationCSR_SeqAIJ", 3379 MatSeqAIJSetPreallocationCSR_SeqAIJ);CHKERRQ(ierr); 3380 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatReorderForNonzeroDiagonal_C", 3381 "MatReorderForNonzeroDiagonal_SeqAIJ", 3382 MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 3383 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMult_seqdense_seqaij_C", 3384 "MatMatMult_SeqDense_SeqAIJ", 3385 MatMatMult_SeqDense_SeqAIJ);CHKERRQ(ierr); 3386 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultSymbolic_seqdense_seqaij_C", 3387 "MatMatMultSymbolic_SeqDense_SeqAIJ", 3388 MatMatMultSymbolic_SeqDense_SeqAIJ);CHKERRQ(ierr); 3389 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultNumeric_seqdense_seqaij_C", 3390 "MatMatMultNumeric_SeqDense_SeqAIJ", 3391 MatMatMultNumeric_SeqDense_SeqAIJ);CHKERRQ(ierr); 3392 ierr = MatCreate_SeqAIJ_Inode(B);CHKERRQ(ierr); 3393 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3394 PetscFunctionReturn(0); 3395 } 3396 EXTERN_C_END 3397 3398 #undef __FUNCT__ 3399 #define __FUNCT__ "MatDuplicateNoCreate_SeqAIJ" 3400 /* 3401 Given a matrix generated with MatGetFactor() duplicates all the information in A into B 3402 */ 3403 PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscTruth mallocmatspace) 3404 { 3405 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 3406 PetscErrorCode ierr; 3407 PetscInt i,m = A->rmap->n; 3408 3409 PetscFunctionBegin; 3410 c = (Mat_SeqAIJ*)C->data; 3411 3412 C->factor = A->factor; 3413 3414 c->row = 0; 3415 c->col = 0; 3416 c->icol = 0; 3417 c->reallocs = 0; 3418 3419 C->assembled = PETSC_TRUE; 3420 3421 ierr = PetscLayoutSetBlockSize(C->rmap,1);CHKERRQ(ierr); 3422 ierr = PetscLayoutSetBlockSize(C->cmap,1);CHKERRQ(ierr); 3423 ierr = PetscLayoutSetUp(C->rmap);CHKERRQ(ierr); 3424 ierr = PetscLayoutSetUp(C->cmap);CHKERRQ(ierr); 3425 3426 ierr = PetscMalloc2(m,PetscInt,&c->imax,m,PetscInt,&c->ilen);CHKERRQ(ierr); 3427 ierr = PetscLogObjectMemory(C, 2*m*sizeof(PetscInt));CHKERRQ(ierr); 3428 for (i=0; i<m; i++) { 3429 c->imax[i] = a->imax[i]; 3430 c->ilen[i] = a->ilen[i]; 3431 } 3432 3433 /* allocate the matrix space */ 3434 if (mallocmatspace){ 3435 ierr = PetscMalloc3(a->i[m],PetscScalar,&c->a,a->i[m],PetscInt,&c->j,m+1,PetscInt,&c->i);CHKERRQ(ierr); 3436 ierr = PetscLogObjectMemory(C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3437 c->singlemalloc = PETSC_TRUE; 3438 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3439 if (m > 0) { 3440 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(PetscInt));CHKERRQ(ierr); 3441 if (cpvalues == MAT_COPY_VALUES) { 3442 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 3443 } else { 3444 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 3445 } 3446 } 3447 } 3448 3449 c->ignorezeroentries = a->ignorezeroentries; 3450 c->roworiented = a->roworiented; 3451 c->nonew = a->nonew; 3452 if (a->diag) { 3453 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 3454 ierr = PetscLogObjectMemory(C,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3455 for (i=0; i<m; i++) { 3456 c->diag[i] = a->diag[i]; 3457 } 3458 } else c->diag = 0; 3459 c->solve_work = 0; 3460 c->saved_values = 0; 3461 c->idiag = 0; 3462 c->ssor_work = 0; 3463 c->keepnonzeropattern = a->keepnonzeropattern; 3464 c->free_a = PETSC_TRUE; 3465 c->free_ij = PETSC_TRUE; 3466 c->xtoy = 0; 3467 c->XtoY = 0; 3468 3469 c->nz = a->nz; 3470 c->maxnz = a->maxnz; 3471 C->preallocated = PETSC_TRUE; 3472 3473 c->compressedrow.use = a->compressedrow.use; 3474 c->compressedrow.nrows = a->compressedrow.nrows; 3475 c->compressedrow.checked = a->compressedrow.checked; 3476 if (a->compressedrow.checked && a->compressedrow.use){ 3477 i = a->compressedrow.nrows; 3478 ierr = PetscMalloc2(i+1,PetscInt,&c->compressedrow.i,i,PetscInt,&c->compressedrow.rindex);CHKERRQ(ierr); 3479 ierr = PetscMemcpy(c->compressedrow.i,a->compressedrow.i,(i+1)*sizeof(PetscInt));CHKERRQ(ierr); 3480 ierr = PetscMemcpy(c->compressedrow.rindex,a->compressedrow.rindex,i*sizeof(PetscInt));CHKERRQ(ierr); 3481 } else { 3482 c->compressedrow.use = PETSC_FALSE; 3483 c->compressedrow.i = PETSC_NULL; 3484 c->compressedrow.rindex = PETSC_NULL; 3485 } 3486 C->same_nonzero = A->same_nonzero; 3487 ierr = MatDuplicate_SeqAIJ_Inode(A,cpvalues,&C);CHKERRQ(ierr); 3488 3489 ierr = PetscFListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 3490 PetscFunctionReturn(0); 3491 } 3492 3493 #undef __FUNCT__ 3494 #define __FUNCT__ "MatDuplicate_SeqAIJ" 3495 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 3496 { 3497 PetscErrorCode ierr; 3498 3499 PetscFunctionBegin; 3500 ierr = MatCreate(((PetscObject)A)->comm,B);CHKERRQ(ierr); 3501 ierr = MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 3502 ierr = MatSetType(*B,MATSEQAIJ);CHKERRQ(ierr); 3503 ierr = MatDuplicateNoCreate_SeqAIJ(*B,A,cpvalues,PETSC_TRUE);CHKERRQ(ierr); 3504 PetscFunctionReturn(0); 3505 } 3506 3507 #undef __FUNCT__ 3508 #define __FUNCT__ "MatLoad_SeqAIJ" 3509 PetscErrorCode MatLoad_SeqAIJ(PetscViewer viewer, const MatType type,Mat *A) 3510 { 3511 Mat_SeqAIJ *a; 3512 Mat B; 3513 PetscErrorCode ierr; 3514 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N; 3515 int fd; 3516 PetscMPIInt size; 3517 MPI_Comm comm; 3518 3519 PetscFunctionBegin; 3520 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 3521 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 3522 if (size > 1) SETERRQ(PETSC_ERR_ARG_SIZ,"view must have one processor"); 3523 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 3524 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 3525 if (header[0] != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 3526 M = header[1]; N = header[2]; nz = header[3]; 3527 3528 if (nz < 0) { 3529 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 3530 } 3531 3532 /* read in row lengths */ 3533 ierr = PetscMalloc(M*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 3534 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 3535 3536 /* check if sum of rowlengths is same as nz */ 3537 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 3538 if (sum != nz) SETERRQ2(PETSC_ERR_FILE_READ,"Inconsistant matrix data in file. no-nonzeros = %d, sum-row-lengths = %d\n",nz,sum); 3539 3540 /* create our matrix */ 3541 ierr = MatCreate(comm,&B);CHKERRQ(ierr); 3542 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 3543 ierr = MatSetType(B,type);CHKERRQ(ierr); 3544 ierr = MatSeqAIJSetPreallocation_SeqAIJ(B,0,rowlengths);CHKERRQ(ierr); 3545 a = (Mat_SeqAIJ*)B->data; 3546 3547 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 3548 3549 /* read in nonzero values */ 3550 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 3551 3552 /* set matrix "i" values */ 3553 a->i[0] = 0; 3554 for (i=1; i<= M; i++) { 3555 a->i[i] = a->i[i-1] + rowlengths[i-1]; 3556 a->ilen[i-1] = rowlengths[i-1]; 3557 } 3558 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 3559 3560 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3561 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3562 *A = B; 3563 PetscFunctionReturn(0); 3564 } 3565 3566 #undef __FUNCT__ 3567 #define __FUNCT__ "MatEqual_SeqAIJ" 3568 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg) 3569 { 3570 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data,*b = (Mat_SeqAIJ *)B->data; 3571 PetscErrorCode ierr; 3572 #if defined(PETSC_USE_COMPLEX) 3573 PetscInt k; 3574 #endif 3575 3576 PetscFunctionBegin; 3577 /* If the matrix dimensions are not equal,or no of nonzeros */ 3578 if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) ||(a->nz != b->nz)) { 3579 *flg = PETSC_FALSE; 3580 PetscFunctionReturn(0); 3581 } 3582 3583 /* if the a->i are the same */ 3584 ierr = PetscMemcmp(a->i,b->i,(A->rmap->n+1)*sizeof(PetscInt),flg);CHKERRQ(ierr); 3585 if (!*flg) PetscFunctionReturn(0); 3586 3587 /* if a->j are the same */ 3588 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),flg);CHKERRQ(ierr); 3589 if (!*flg) PetscFunctionReturn(0); 3590 3591 /* if a->a are the same */ 3592 #if defined(PETSC_USE_COMPLEX) 3593 for (k=0; k<a->nz; k++){ 3594 if (PetscRealPart(a->a[k]) != PetscRealPart(b->a[k]) || PetscImaginaryPart(a->a[k]) != PetscImaginaryPart(b->a[k])){ 3595 *flg = PETSC_FALSE; 3596 PetscFunctionReturn(0); 3597 } 3598 } 3599 #else 3600 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 3601 #endif 3602 PetscFunctionReturn(0); 3603 } 3604 3605 #undef __FUNCT__ 3606 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 3607 /*@ 3608 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 3609 provided by the user. 3610 3611 Collective on MPI_Comm 3612 3613 Input Parameters: 3614 + comm - must be an MPI communicator of size 1 3615 . m - number of rows 3616 . n - number of columns 3617 . i - row indices 3618 . j - column indices 3619 - a - matrix values 3620 3621 Output Parameter: 3622 . mat - the matrix 3623 3624 Level: intermediate 3625 3626 Notes: 3627 The i, j, and a arrays are not copied by this routine, the user must free these arrays 3628 once the matrix is destroyed 3629 3630 You cannot set new nonzero locations into this matrix, that will generate an error. 3631 3632 The i and j indices are 0 based 3633 3634 The format which is used for the sparse matrix input, is equivalent to a 3635 row-major ordering.. i.e for the following matrix, the input data expected is 3636 as shown: 3637 3638 1 0 0 3639 2 0 3 3640 4 5 6 3641 3642 i = {0,1,3,6} [size = nrow+1 = 3+1] 3643 j = {0,0,2,0,1,2} [size = nz = 6]; values must be sorted for each row 3644 v = {1,2,3,4,5,6} [size = nz = 6] 3645 3646 3647 .seealso: MatCreate(), MatCreateMPIAIJ(), MatCreateSeqAIJ(), MatCreateMPIAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 3648 3649 @*/ 3650 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt* i,PetscInt*j,PetscScalar *a,Mat *mat) 3651 { 3652 PetscErrorCode ierr; 3653 PetscInt ii; 3654 Mat_SeqAIJ *aij; 3655 #if defined(PETSC_USE_DEBUG) 3656 PetscInt jj; 3657 #endif 3658 3659 PetscFunctionBegin; 3660 if (i[0]) { 3661 SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 3662 } 3663 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 3664 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 3665 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 3666 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 3667 aij = (Mat_SeqAIJ*)(*mat)->data; 3668 ierr = PetscMalloc2(m,PetscInt,&aij->imax,m,PetscInt,&aij->ilen);CHKERRQ(ierr); 3669 3670 aij->i = i; 3671 aij->j = j; 3672 aij->a = a; 3673 aij->singlemalloc = PETSC_FALSE; 3674 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 3675 aij->free_a = PETSC_FALSE; 3676 aij->free_ij = PETSC_FALSE; 3677 3678 for (ii=0; ii<m; ii++) { 3679 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 3680 #if defined(PETSC_USE_DEBUG) 3681 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]); 3682 for (jj=i[ii]+1; jj<i[ii+1]; jj++) { 3683 if (j[jj] < j[jj-1]) SETERRQ3(PETSC_ERR_ARG_OUTOFRANGE,"Column entry number %D (actual colum %D) in row %D is not sorted",jj-i[ii],j[jj],ii); 3684 if (j[jj] == j[jj]-1) SETERRQ3(PETSC_ERR_ARG_OUTOFRANGE,"Column entry number %D (actual colum %D) in row %D is identical to previous entry",jj-i[ii],j[jj],ii); 3685 } 3686 #endif 3687 } 3688 #if defined(PETSC_USE_DEBUG) 3689 for (ii=0; ii<aij->i[m]; ii++) { 3690 if (j[ii] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 3691 if (j[ii] > n - 1) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column index to large at location = %d index = %d",ii,j[ii]); 3692 } 3693 #endif 3694 3695 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3696 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3697 PetscFunctionReturn(0); 3698 } 3699 3700 #undef __FUNCT__ 3701 #define __FUNCT__ "MatSetColoring_SeqAIJ" 3702 PetscErrorCode MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 3703 { 3704 PetscErrorCode ierr; 3705 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3706 3707 PetscFunctionBegin; 3708 if (coloring->ctype == IS_COLORING_GLOBAL) { 3709 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 3710 a->coloring = coloring; 3711 } else if (coloring->ctype == IS_COLORING_GHOSTED) { 3712 PetscInt i,*larray; 3713 ISColoring ocoloring; 3714 ISColoringValue *colors; 3715 3716 /* set coloring for diagonal portion */ 3717 ierr = PetscMalloc(A->cmap->n*sizeof(PetscInt),&larray);CHKERRQ(ierr); 3718 for (i=0; i<A->cmap->n; i++) { 3719 larray[i] = i; 3720 } 3721 ierr = ISGlobalToLocalMappingApply(A->mapping,IS_GTOLM_MASK,A->cmap->n,larray,PETSC_NULL,larray);CHKERRQ(ierr); 3722 ierr = PetscMalloc(A->cmap->n*sizeof(ISColoringValue),&colors);CHKERRQ(ierr); 3723 for (i=0; i<A->cmap->n; i++) { 3724 colors[i] = coloring->colors[larray[i]]; 3725 } 3726 ierr = PetscFree(larray);CHKERRQ(ierr); 3727 ierr = ISColoringCreate(PETSC_COMM_SELF,coloring->n,A->cmap->n,colors,&ocoloring);CHKERRQ(ierr); 3728 a->coloring = ocoloring; 3729 } 3730 PetscFunctionReturn(0); 3731 } 3732 3733 #if defined(PETSC_HAVE_ADIC) 3734 EXTERN_C_BEGIN 3735 #include "adic/ad_utils.h" 3736 EXTERN_C_END 3737 3738 #undef __FUNCT__ 3739 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 3740 PetscErrorCode MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 3741 { 3742 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3743 PetscInt m = A->rmap->n,*ii = a->i,*jj = a->j,nz,i,j,nlen; 3744 PetscScalar *v = a->a,*values = ((PetscScalar*)advalues)+1; 3745 ISColoringValue *color; 3746 3747 PetscFunctionBegin; 3748 if (!a->coloring) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 3749 nlen = PetscADGetDerivTypeSize()/sizeof(PetscScalar); 3750 color = a->coloring->colors; 3751 /* loop over rows */ 3752 for (i=0; i<m; i++) { 3753 nz = ii[i+1] - ii[i]; 3754 /* loop over columns putting computed value into matrix */ 3755 for (j=0; j<nz; j++) { 3756 *v++ = values[color[*jj++]]; 3757 } 3758 values += nlen; /* jump to next row of derivatives */ 3759 } 3760 PetscFunctionReturn(0); 3761 } 3762 #endif 3763 3764 #undef __FUNCT__ 3765 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 3766 PetscErrorCode MatSetValuesAdifor_SeqAIJ(Mat A,PetscInt nl,void *advalues) 3767 { 3768 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3769 PetscInt m = A->rmap->n,*ii = a->i,*jj = a->j,nz,i,j; 3770 MatScalar *v = a->a; 3771 PetscScalar *values = (PetscScalar *)advalues; 3772 ISColoringValue *color; 3773 3774 PetscFunctionBegin; 3775 if (!a->coloring) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 3776 color = a->coloring->colors; 3777 /* loop over rows */ 3778 for (i=0; i<m; i++) { 3779 nz = ii[i+1] - ii[i]; 3780 /* loop over columns putting computed value into matrix */ 3781 for (j=0; j<nz; j++) { 3782 *v++ = values[color[*jj++]]; 3783 } 3784 values += nl; /* jump to next row of derivatives */ 3785 } 3786 PetscFunctionReturn(0); 3787 } 3788 3789 /* 3790 Special version for direct calls from Fortran 3791 */ 3792 #if defined(PETSC_HAVE_FORTRAN_CAPS) 3793 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 3794 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 3795 #define matsetvaluesseqaij_ matsetvaluesseqaij 3796 #endif 3797 3798 /* Change these macros so can be used in void function */ 3799 #undef CHKERRQ 3800 #define CHKERRQ(ierr) CHKERRABORT(((PetscObject)A)->comm,ierr) 3801 #undef SETERRQ2 3802 #define SETERRQ2(ierr,b,c,d) CHKERRABORT(((PetscObject)A)->comm,ierr) 3803 3804 EXTERN_C_BEGIN 3805 #undef __FUNCT__ 3806 #define __FUNCT__ "matsetvaluesseqaij_" 3807 void PETSC_STDCALL matsetvaluesseqaij_(Mat *AA,PetscInt *mm,const PetscInt im[],PetscInt *nn,const PetscInt in[],const PetscScalar v[],InsertMode *isis, PetscErrorCode *_ierr) 3808 { 3809 Mat A = *AA; 3810 PetscInt m = *mm, n = *nn; 3811 InsertMode is = *isis; 3812 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3813 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 3814 PetscInt *imax,*ai,*ailen; 3815 PetscErrorCode ierr; 3816 PetscInt *aj,nonew = a->nonew,lastcol = -1; 3817 MatScalar *ap,value,*aa; 3818 PetscTruth ignorezeroentries = a->ignorezeroentries; 3819 PetscTruth roworiented = a->roworiented; 3820 3821 PetscFunctionBegin; 3822 ierr = MatPreallocated(A);CHKERRQ(ierr); 3823 imax = a->imax; 3824 ai = a->i; 3825 ailen = a->ilen; 3826 aj = a->j; 3827 aa = a->a; 3828 3829 for (k=0; k<m; k++) { /* loop over added rows */ 3830 row = im[k]; 3831 if (row < 0) continue; 3832 #if defined(PETSC_USE_DEBUG) 3833 if (row >= A->rmap->n) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 3834 #endif 3835 rp = aj + ai[row]; ap = aa + ai[row]; 3836 rmax = imax[row]; nrow = ailen[row]; 3837 low = 0; 3838 high = nrow; 3839 for (l=0; l<n; l++) { /* loop over added columns */ 3840 if (in[l] < 0) continue; 3841 #if defined(PETSC_USE_DEBUG) 3842 if (in[l] >= A->cmap->n) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 3843 #endif 3844 col = in[l]; 3845 if (roworiented) { 3846 value = v[l + k*n]; 3847 } else { 3848 value = v[k + l*m]; 3849 } 3850 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 3851 3852 if (col <= lastcol) low = 0; else high = nrow; 3853 lastcol = col; 3854 while (high-low > 5) { 3855 t = (low+high)/2; 3856 if (rp[t] > col) high = t; 3857 else low = t; 3858 } 3859 for (i=low; i<high; i++) { 3860 if (rp[i] > col) break; 3861 if (rp[i] == col) { 3862 if (is == ADD_VALUES) ap[i] += value; 3863 else ap[i] = value; 3864 goto noinsert; 3865 } 3866 } 3867 if (value == 0.0 && ignorezeroentries) goto noinsert; 3868 if (nonew == 1) goto noinsert; 3869 if (nonew == -1) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 3870 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 3871 N = nrow++ - 1; a->nz++; high++; 3872 /* shift up all the later entries in this row */ 3873 for (ii=N; ii>=i; ii--) { 3874 rp[ii+1] = rp[ii]; 3875 ap[ii+1] = ap[ii]; 3876 } 3877 rp[i] = col; 3878 ap[i] = value; 3879 noinsert:; 3880 low = i + 1; 3881 } 3882 ailen[row] = nrow; 3883 } 3884 A->same_nonzero = PETSC_FALSE; 3885 PetscFunctionReturnVoid(); 3886 } 3887 EXTERN_C_END 3888