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