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