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