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.use){ierr = PetscFree(a->compressedrow.i);} 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 = PetscMalloc((2*nrows+1)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 1695 starts = lens + nrows; 1696 /* loop over new rows determining lens and starting points */ 1697 for (i=0; i<nrows; i++) { 1698 kstart = ai[irow[i]]; 1699 kend = kstart + ailen[irow[i]]; 1700 for (k=kstart; k<kend; k++) { 1701 if (aj[k] >= first) { 1702 starts[i] = k; 1703 break; 1704 } 1705 } 1706 sum = 0; 1707 while (k < kend) { 1708 if (aj[k++] >= first+ncols) break; 1709 sum++; 1710 } 1711 lens[i] = sum; 1712 } 1713 /* create submatrix */ 1714 if (scall == MAT_REUSE_MATRIX) { 1715 PetscInt n_cols,n_rows; 1716 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 1717 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 1718 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 1719 C = *B; 1720 } else { 1721 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 1722 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1723 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 1724 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 1725 } 1726 c = (Mat_SeqAIJ*)C->data; 1727 1728 /* loop over rows inserting into submatrix */ 1729 a_new = c->a; 1730 j_new = c->j; 1731 i_new = c->i; 1732 1733 for (i=0; i<nrows; i++) { 1734 ii = starts[i]; 1735 lensi = lens[i]; 1736 for (k=0; k<lensi; k++) { 1737 *j_new++ = aj[ii+k] - first; 1738 } 1739 ierr = PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(PetscScalar));CHKERRQ(ierr); 1740 a_new += lensi; 1741 i_new[i+1] = i_new[i] + lensi; 1742 c->ilen[i] = lensi; 1743 } 1744 ierr = PetscFree(lens);CHKERRQ(ierr); 1745 } else { 1746 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 1747 ierr = PetscMalloc((1+oldcols)*sizeof(PetscInt),&smap);CHKERRQ(ierr); 1748 1749 ierr = PetscMalloc((1+nrows)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 1750 ierr = PetscMemzero(smap,oldcols*sizeof(PetscInt));CHKERRQ(ierr); 1751 for (i=0; i<ncols; i++) smap[icol[i]] = i+1; 1752 /* determine lens of each row */ 1753 for (i=0; i<nrows; i++) { 1754 kstart = ai[irow[i]]; 1755 kend = kstart + a->ilen[irow[i]]; 1756 lens[i] = 0; 1757 for (k=kstart; k<kend; k++) { 1758 if (smap[aj[k]]) { 1759 lens[i]++; 1760 } 1761 } 1762 } 1763 /* Create and fill new matrix */ 1764 if (scall == MAT_REUSE_MATRIX) { 1765 PetscTruth equal; 1766 1767 c = (Mat_SeqAIJ *)((*B)->data); 1768 if ((*B)->rmap->n != nrows || (*B)->cmap->n != ncols) SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 1769 ierr = PetscMemcmp(c->ilen,lens,(*B)->rmap->n*sizeof(PetscInt),&equal);CHKERRQ(ierr); 1770 if (!equal) { 1771 SETERRQ(PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 1772 } 1773 ierr = PetscMemzero(c->ilen,(*B)->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 1774 C = *B; 1775 } else { 1776 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 1777 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 1778 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 1779 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 1780 } 1781 c = (Mat_SeqAIJ *)(C->data); 1782 for (i=0; i<nrows; i++) { 1783 row = irow[i]; 1784 kstart = ai[row]; 1785 kend = kstart + a->ilen[row]; 1786 mat_i = c->i[i]; 1787 mat_j = c->j + mat_i; 1788 mat_a = c->a + mat_i; 1789 mat_ilen = c->ilen + i; 1790 for (k=kstart; k<kend; k++) { 1791 if ((tcol=smap[a->j[k]])) { 1792 *mat_j++ = tcol - 1; 1793 *mat_a++ = a->a[k]; 1794 (*mat_ilen)++; 1795 1796 } 1797 } 1798 } 1799 /* Free work space */ 1800 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 1801 ierr = PetscFree(smap);CHKERRQ(ierr); 1802 ierr = PetscFree(lens);CHKERRQ(ierr); 1803 } 1804 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1805 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1806 1807 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 1808 *B = C; 1809 PetscFunctionReturn(0); 1810 } 1811 1812 #undef __FUNCT__ 1813 #define __FUNCT__ "MatILUFactor_SeqAIJ_newdatastruct" 1814 PetscErrorCode MatILUFactor_SeqAIJ_newdatastruct(Mat inA,IS row,IS col,const MatFactorInfo *info) 1815 { 1816 PetscErrorCode ierr; 1817 Mat C; 1818 1819 PetscFunctionBegin; 1820 /* A new factored matrix C replaces original matrix inA, not a truely in-place ILUFactor yet */ 1821 ierr = MatGetFactor(inA,MAT_SOLVER_PETSC,MAT_FACTOR_LU,&C);CHKERRQ(ierr); 1822 ierr = MatILUFactorSymbolic_SeqAIJ_newdatastruct(C,inA,row,col,info);CHKERRQ(ierr); 1823 ierr = MatLUFactorNumeric_SeqAIJ_newdatastruct(C,inA,info);CHKERRQ(ierr); 1824 inA->ops->solve = C->ops->solve; 1825 inA->ops->solvetranspose = C->ops->solvetranspose; 1826 ierr = MatHeaderCopy(inA,C);CHKERRQ(ierr); 1827 ierr = PetscLogObjectParent(inA,((Mat_SeqAIJ*)(inA->data))->icol);CHKERRQ(ierr); 1828 PetscFunctionReturn(0); 1829 } 1830 1831 #undef __FUNCT__ 1832 #define __FUNCT__ "MatILUFactor_SeqAIJ" 1833 PetscErrorCode MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,const MatFactorInfo *info) 1834 { 1835 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1836 PetscErrorCode ierr; 1837 Mat outA; 1838 PetscTruth row_identity,col_identity; 1839 PetscTruth newdatastruct=PETSC_FALSE; 1840 1841 PetscFunctionBegin; 1842 ierr = PetscOptionsGetTruth(PETSC_NULL,"-ilu_new",&newdatastruct,PETSC_NULL);CHKERRQ(ierr); 1843 if(newdatastruct){ 1844 ierr = MatILUFactor_SeqAIJ_newdatastruct(inA,row,col,info);CHKERRQ(ierr); 1845 PetscFunctionReturn(0); 1846 } 1847 1848 if (info->levels != 0) SETERRQ(PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 1849 1850 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 1851 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 1852 1853 outA = inA; 1854 outA->factor = MAT_FACTOR_LU; 1855 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1856 if (a->row) { ierr = ISDestroy(a->row);CHKERRQ(ierr);} 1857 a->row = row; 1858 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 1859 if (a->col) { ierr = ISDestroy(a->col);CHKERRQ(ierr);} 1860 a->col = col; 1861 1862 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 1863 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} /* need to remove old one */ 1864 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 1865 ierr = PetscLogObjectParent(inA,a->icol);CHKERRQ(ierr); 1866 1867 if (!a->solve_work) { /* this matrix may have been factored before */ 1868 ierr = PetscMalloc((inA->rmap->n+1)*sizeof(PetscScalar),&a->solve_work);CHKERRQ(ierr); 1869 ierr = PetscLogObjectMemory(inA, (inA->rmap->n+1)*sizeof(PetscScalar));CHKERRQ(ierr); 1870 } 1871 1872 ierr = MatMarkDiagonal_SeqAIJ(inA);CHKERRQ(ierr); 1873 if (row_identity && col_identity) { 1874 ierr = MatLUFactorNumeric_SeqAIJ(outA,inA,info);CHKERRQ(ierr); 1875 } else { 1876 ierr = MatLUFactorNumeric_SeqAIJ_InplaceWithPerm(outA,inA,info);CHKERRQ(ierr); 1877 } 1878 PetscFunctionReturn(0); 1879 } 1880 1881 #include "petscblaslapack.h" 1882 #undef __FUNCT__ 1883 #define __FUNCT__ "MatScale_SeqAIJ" 1884 PetscErrorCode MatScale_SeqAIJ(Mat inA,PetscScalar alpha) 1885 { 1886 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 1887 PetscScalar oalpha = alpha; 1888 PetscErrorCode ierr; 1889 PetscBLASInt one = 1,bnz = PetscBLASIntCast(a->nz); 1890 1891 PetscFunctionBegin; 1892 BLASscal_(&bnz,&oalpha,a->a,&one); 1893 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1894 PetscFunctionReturn(0); 1895 } 1896 1897 #undef __FUNCT__ 1898 #define __FUNCT__ "MatGetSubMatrices_SeqAIJ" 1899 PetscErrorCode MatGetSubMatrices_SeqAIJ(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 1900 { 1901 PetscErrorCode ierr; 1902 PetscInt i; 1903 1904 PetscFunctionBegin; 1905 if (scall == MAT_INITIAL_MATRIX) { 1906 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 1907 } 1908 1909 for (i=0; i<n; i++) { 1910 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 1911 } 1912 PetscFunctionReturn(0); 1913 } 1914 1915 #undef __FUNCT__ 1916 #define __FUNCT__ "MatIncreaseOverlap_SeqAIJ" 1917 PetscErrorCode MatIncreaseOverlap_SeqAIJ(Mat A,PetscInt is_max,IS is[],PetscInt ov) 1918 { 1919 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1920 PetscErrorCode ierr; 1921 PetscInt row,i,j,k,l,m,n,*nidx,isz,val; 1922 const PetscInt *idx; 1923 PetscInt start,end,*ai,*aj; 1924 PetscBT table; 1925 1926 PetscFunctionBegin; 1927 m = A->rmap->n; 1928 ai = a->i; 1929 aj = a->j; 1930 1931 if (ov < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"illegal negative overlap value used"); 1932 1933 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&nidx);CHKERRQ(ierr); 1934 ierr = PetscBTCreate(m,table);CHKERRQ(ierr); 1935 1936 for (i=0; i<is_max; i++) { 1937 /* Initialize the two local arrays */ 1938 isz = 0; 1939 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 1940 1941 /* Extract the indices, assume there can be duplicate entries */ 1942 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 1943 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 1944 1945 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 1946 for (j=0; j<n ; ++j){ 1947 if(!PetscBTLookupSet(table,idx[j])) { nidx[isz++] = idx[j];} 1948 } 1949 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 1950 ierr = ISDestroy(is[i]);CHKERRQ(ierr); 1951 1952 k = 0; 1953 for (j=0; j<ov; j++){ /* for each overlap */ 1954 n = isz; 1955 for (; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 1956 row = nidx[k]; 1957 start = ai[row]; 1958 end = ai[row+1]; 1959 for (l = start; l<end ; l++){ 1960 val = aj[l] ; 1961 if (!PetscBTLookupSet(table,val)) {nidx[isz++] = val;} 1962 } 1963 } 1964 } 1965 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,(is+i));CHKERRQ(ierr); 1966 } 1967 ierr = PetscBTDestroy(table);CHKERRQ(ierr); 1968 ierr = PetscFree(nidx);CHKERRQ(ierr); 1969 PetscFunctionReturn(0); 1970 } 1971 1972 /* -------------------------------------------------------------- */ 1973 #undef __FUNCT__ 1974 #define __FUNCT__ "MatPermute_SeqAIJ" 1975 PetscErrorCode MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 1976 { 1977 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1978 PetscErrorCode ierr; 1979 PetscInt i,nz,m = A->rmap->n,n = A->cmap->n; 1980 const PetscInt *row,*col; 1981 PetscInt *cnew,j,*lens; 1982 IS icolp,irowp; 1983 PetscInt *cwork; 1984 PetscScalar *vwork; 1985 1986 PetscFunctionBegin; 1987 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 1988 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 1989 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 1990 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 1991 1992 /* determine lengths of permuted rows */ 1993 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 1994 for (i=0; i<m; i++) { 1995 lens[row[i]] = a->i[i+1] - a->i[i]; 1996 } 1997 ierr = MatCreate(((PetscObject)A)->comm,B);CHKERRQ(ierr); 1998 ierr = MatSetSizes(*B,m,n,m,n);CHKERRQ(ierr); 1999 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 2000 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*B,0,lens);CHKERRQ(ierr); 2001 ierr = PetscFree(lens);CHKERRQ(ierr); 2002 2003 ierr = PetscMalloc(n*sizeof(PetscInt),&cnew);CHKERRQ(ierr); 2004 for (i=0; i<m; i++) { 2005 ierr = MatGetRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2006 for (j=0; j<nz; j++) { cnew[j] = col[cwork[j]];} 2007 ierr = MatSetValues_SeqAIJ(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 2008 ierr = MatRestoreRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2009 } 2010 ierr = PetscFree(cnew);CHKERRQ(ierr); 2011 (*B)->assembled = PETSC_FALSE; 2012 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2013 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2014 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 2015 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 2016 ierr = ISDestroy(irowp);CHKERRQ(ierr); 2017 ierr = ISDestroy(icolp);CHKERRQ(ierr); 2018 PetscFunctionReturn(0); 2019 } 2020 2021 #undef __FUNCT__ 2022 #define __FUNCT__ "MatCopy_SeqAIJ" 2023 PetscErrorCode MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 2024 { 2025 PetscErrorCode ierr; 2026 2027 PetscFunctionBegin; 2028 /* If the two matrices have the same copy implementation, use fast copy. */ 2029 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 2030 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2031 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 2032 2033 if (a->i[A->rmap->n] != b->i[B->rmap->n]) { 2034 SETERRQ(PETSC_ERR_ARG_INCOMP,"Number of nonzeros in two matrices are different"); 2035 } 2036 ierr = PetscMemcpy(b->a,a->a,(a->i[A->rmap->n])*sizeof(PetscScalar));CHKERRQ(ierr); 2037 } else { 2038 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 2039 } 2040 PetscFunctionReturn(0); 2041 } 2042 2043 #undef __FUNCT__ 2044 #define __FUNCT__ "MatSetUpPreallocation_SeqAIJ" 2045 PetscErrorCode MatSetUpPreallocation_SeqAIJ(Mat A) 2046 { 2047 PetscErrorCode ierr; 2048 2049 PetscFunctionBegin; 2050 ierr = MatSeqAIJSetPreallocation_SeqAIJ(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 2051 PetscFunctionReturn(0); 2052 } 2053 2054 #undef __FUNCT__ 2055 #define __FUNCT__ "MatGetArray_SeqAIJ" 2056 PetscErrorCode MatGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 2057 { 2058 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2059 PetscFunctionBegin; 2060 *array = a->a; 2061 PetscFunctionReturn(0); 2062 } 2063 2064 #undef __FUNCT__ 2065 #define __FUNCT__ "MatRestoreArray_SeqAIJ" 2066 PetscErrorCode MatRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 2067 { 2068 PetscFunctionBegin; 2069 PetscFunctionReturn(0); 2070 } 2071 2072 #undef __FUNCT__ 2073 #define __FUNCT__ "MatFDColoringApply_SeqAIJ" 2074 PetscErrorCode MatFDColoringApply_SeqAIJ(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx) 2075 { 2076 PetscErrorCode (*f)(void*,Vec,Vec,void*) = (PetscErrorCode (*)(void*,Vec,Vec,void *))coloring->f; 2077 PetscErrorCode ierr; 2078 PetscInt k,N,start,end,l,row,col,srow,**vscaleforrow,m1,m2; 2079 PetscScalar dx,*y,*xx,*w3_array; 2080 PetscScalar *vscale_array; 2081 PetscReal epsilon = coloring->error_rel,umin = coloring->umin; 2082 Vec w1,w2,w3; 2083 void *fctx = coloring->fctx; 2084 PetscTruth flg = PETSC_FALSE; 2085 2086 PetscFunctionBegin; 2087 if (!coloring->w1) { 2088 ierr = VecDuplicate(x1,&coloring->w1);CHKERRQ(ierr); 2089 ierr = PetscLogObjectParent(coloring,coloring->w1);CHKERRQ(ierr); 2090 ierr = VecDuplicate(x1,&coloring->w2);CHKERRQ(ierr); 2091 ierr = PetscLogObjectParent(coloring,coloring->w2);CHKERRQ(ierr); 2092 ierr = VecDuplicate(x1,&coloring->w3);CHKERRQ(ierr); 2093 ierr = PetscLogObjectParent(coloring,coloring->w3);CHKERRQ(ierr); 2094 } 2095 w1 = coloring->w1; w2 = coloring->w2; w3 = coloring->w3; 2096 2097 ierr = MatSetUnfactored(J);CHKERRQ(ierr); 2098 ierr = PetscOptionsGetTruth(((PetscObject)coloring)->prefix,"-mat_fd_coloring_dont_rezero",&flg,PETSC_NULL);CHKERRQ(ierr); 2099 if (flg) { 2100 ierr = PetscInfo(coloring,"Not calling MatZeroEntries()\n");CHKERRQ(ierr); 2101 } else { 2102 PetscTruth assembled; 2103 ierr = MatAssembled(J,&assembled);CHKERRQ(ierr); 2104 if (assembled) { 2105 ierr = MatZeroEntries(J);CHKERRQ(ierr); 2106 } 2107 } 2108 2109 ierr = VecGetOwnershipRange(x1,&start,&end);CHKERRQ(ierr); 2110 ierr = VecGetSize(x1,&N);CHKERRQ(ierr); 2111 2112 /* 2113 This is a horrible, horrible, hack. See DMMGComputeJacobian_Multigrid() it inproperly sets 2114 coloring->F for the coarser grids from the finest 2115 */ 2116 if (coloring->F) { 2117 ierr = VecGetLocalSize(coloring->F,&m1);CHKERRQ(ierr); 2118 ierr = VecGetLocalSize(w1,&m2);CHKERRQ(ierr); 2119 if (m1 != m2) { 2120 coloring->F = 0; 2121 } 2122 } 2123 2124 if (coloring->F) { 2125 w1 = coloring->F; 2126 coloring->F = 0; 2127 } else { 2128 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2129 ierr = (*f)(sctx,x1,w1,fctx);CHKERRQ(ierr); 2130 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2131 } 2132 2133 /* 2134 Compute all the scale factors and share with other processors 2135 */ 2136 ierr = VecGetArray(x1,&xx);CHKERRQ(ierr);xx = xx - start; 2137 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr);vscale_array = vscale_array - start; 2138 for (k=0; k<coloring->ncolors; k++) { 2139 /* 2140 Loop over each column associated with color adding the 2141 perturbation to the vector w3. 2142 */ 2143 for (l=0; l<coloring->ncolumns[k]; l++) { 2144 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2145 dx = xx[col]; 2146 if (dx == 0.0) dx = 1.0; 2147 #if !defined(PETSC_USE_COMPLEX) 2148 if (dx < umin && dx >= 0.0) dx = umin; 2149 else if (dx < 0.0 && dx > -umin) dx = -umin; 2150 #else 2151 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2152 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2153 #endif 2154 dx *= epsilon; 2155 vscale_array[col] = 1.0/dx; 2156 } 2157 } 2158 vscale_array = vscale_array + start;ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2159 ierr = VecGhostUpdateBegin(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2160 ierr = VecGhostUpdateEnd(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2161 2162 /* ierr = VecView(coloring->vscale,PETSC_VIEWER_STDOUT_WORLD); 2163 ierr = VecView(x1,PETSC_VIEWER_STDOUT_WORLD);*/ 2164 2165 if (coloring->vscaleforrow) vscaleforrow = coloring->vscaleforrow; 2166 else vscaleforrow = coloring->columnsforrow; 2167 2168 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2169 /* 2170 Loop over each color 2171 */ 2172 for (k=0; k<coloring->ncolors; k++) { 2173 coloring->currentcolor = k; 2174 ierr = VecCopy(x1,w3);CHKERRQ(ierr); 2175 ierr = VecGetArray(w3,&w3_array);CHKERRQ(ierr);w3_array = w3_array - start; 2176 /* 2177 Loop over each column associated with color adding the 2178 perturbation to the vector w3. 2179 */ 2180 for (l=0; l<coloring->ncolumns[k]; l++) { 2181 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2182 dx = xx[col]; 2183 if (dx == 0.0) dx = 1.0; 2184 #if !defined(PETSC_USE_COMPLEX) 2185 if (dx < umin && dx >= 0.0) dx = umin; 2186 else if (dx < 0.0 && dx > -umin) dx = -umin; 2187 #else 2188 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2189 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2190 #endif 2191 dx *= epsilon; 2192 if (!PetscAbsScalar(dx)) SETERRQ(PETSC_ERR_PLIB,"Computed 0 differencing parameter"); 2193 w3_array[col] += dx; 2194 } 2195 w3_array = w3_array + start; ierr = VecRestoreArray(w3,&w3_array);CHKERRQ(ierr); 2196 2197 /* 2198 Evaluate function at x1 + dx (here dx is a vector of perturbations) 2199 */ 2200 2201 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2202 ierr = (*f)(sctx,w3,w2,fctx);CHKERRQ(ierr); 2203 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2204 ierr = VecAXPY(w2,-1.0,w1);CHKERRQ(ierr); 2205 2206 /* 2207 Loop over rows of vector, putting results into Jacobian matrix 2208 */ 2209 ierr = VecGetArray(w2,&y);CHKERRQ(ierr); 2210 for (l=0; l<coloring->nrows[k]; l++) { 2211 row = coloring->rows[k][l]; 2212 col = coloring->columnsforrow[k][l]; 2213 y[row] *= vscale_array[vscaleforrow[k][l]]; 2214 srow = row + start; 2215 ierr = MatSetValues_SeqAIJ(J,1,&srow,1,&col,y+row,INSERT_VALUES);CHKERRQ(ierr); 2216 } 2217 ierr = VecRestoreArray(w2,&y);CHKERRQ(ierr); 2218 } 2219 coloring->currentcolor = k; 2220 ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2221 xx = xx + start; ierr = VecRestoreArray(x1,&xx);CHKERRQ(ierr); 2222 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2223 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2224 PetscFunctionReturn(0); 2225 } 2226 2227 #include "petscblaslapack.h" 2228 #undef __FUNCT__ 2229 #define __FUNCT__ "MatAXPY_SeqAIJ" 2230 PetscErrorCode MatAXPY_SeqAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 2231 { 2232 PetscErrorCode ierr; 2233 PetscInt i; 2234 Mat_SeqAIJ *x = (Mat_SeqAIJ *)X->data,*y = (Mat_SeqAIJ *)Y->data; 2235 PetscBLASInt one=1,bnz = PetscBLASIntCast(x->nz); 2236 2237 PetscFunctionBegin; 2238 if (str == SAME_NONZERO_PATTERN) { 2239 PetscScalar alpha = a; 2240 BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one); 2241 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 2242 if (y->xtoy && y->XtoY != X) { 2243 ierr = PetscFree(y->xtoy);CHKERRQ(ierr); 2244 ierr = MatDestroy(y->XtoY);CHKERRQ(ierr); 2245 } 2246 if (!y->xtoy) { /* get xtoy */ 2247 ierr = MatAXPYGetxtoy_Private(X->rmap->n,x->i,x->j,PETSC_NULL, y->i,y->j,PETSC_NULL, &y->xtoy);CHKERRQ(ierr); 2248 y->XtoY = X; 2249 ierr = PetscObjectReference((PetscObject)X);CHKERRQ(ierr); 2250 } 2251 for (i=0; i<x->nz; i++) y->a[y->xtoy[i]] += a*(x->a[i]); 2252 ierr = PetscInfo3(Y,"ratio of nnz(X)/nnz(Y): %d/%d = %G\n",x->nz,y->nz,(PetscReal)(x->nz)/y->nz);CHKERRQ(ierr); 2253 } else { 2254 ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr); 2255 } 2256 PetscFunctionReturn(0); 2257 } 2258 2259 #undef __FUNCT__ 2260 #define __FUNCT__ "MatSetBlockSize_SeqAIJ" 2261 PetscErrorCode MatSetBlockSize_SeqAIJ(Mat A,PetscInt bs) 2262 { 2263 PetscErrorCode ierr; 2264 2265 PetscFunctionBegin; 2266 ierr = PetscLayoutSetBlockSize(A->rmap,bs);CHKERRQ(ierr); 2267 ierr = PetscLayoutSetBlockSize(A->cmap,bs);CHKERRQ(ierr); 2268 PetscFunctionReturn(0); 2269 } 2270 2271 #undef __FUNCT__ 2272 #define __FUNCT__ "MatConjugate_SeqAIJ" 2273 PetscErrorCode PETSCMAT_DLLEXPORT MatConjugate_SeqAIJ(Mat mat) 2274 { 2275 #if defined(PETSC_USE_COMPLEX) 2276 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2277 PetscInt i,nz; 2278 PetscScalar *a; 2279 2280 PetscFunctionBegin; 2281 nz = aij->nz; 2282 a = aij->a; 2283 for (i=0; i<nz; i++) { 2284 a[i] = PetscConj(a[i]); 2285 } 2286 #else 2287 PetscFunctionBegin; 2288 #endif 2289 PetscFunctionReturn(0); 2290 } 2291 2292 #undef __FUNCT__ 2293 #define __FUNCT__ "MatGetRowMaxAbs_SeqAIJ" 2294 PetscErrorCode MatGetRowMaxAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2295 { 2296 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2297 PetscErrorCode ierr; 2298 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2299 PetscReal atmp; 2300 PetscScalar *x; 2301 MatScalar *aa; 2302 2303 PetscFunctionBegin; 2304 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2305 aa = a->a; 2306 ai = a->i; 2307 aj = a->j; 2308 2309 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2310 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2311 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2312 if (n != A->rmap->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2313 for (i=0; i<m; i++) { 2314 ncols = ai[1] - ai[0]; ai++; 2315 x[i] = 0.0; 2316 for (j=0; j<ncols; j++){ 2317 atmp = PetscAbsScalar(*aa); 2318 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2319 aa++; aj++; 2320 } 2321 } 2322 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2323 PetscFunctionReturn(0); 2324 } 2325 2326 #undef __FUNCT__ 2327 #define __FUNCT__ "MatGetRowMax_SeqAIJ" 2328 PetscErrorCode MatGetRowMax_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2329 { 2330 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2331 PetscErrorCode ierr; 2332 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2333 PetscScalar *x; 2334 MatScalar *aa; 2335 2336 PetscFunctionBegin; 2337 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2338 aa = a->a; 2339 ai = a->i; 2340 aj = a->j; 2341 2342 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2343 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2344 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2345 if (n != A->rmap->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2346 for (i=0; i<m; i++) { 2347 ncols = ai[1] - ai[0]; ai++; 2348 if (ncols == A->cmap->n) { /* row is dense */ 2349 x[i] = *aa; if (idx) idx[i] = 0; 2350 } else { /* row is sparse so already KNOW maximum is 0.0 or higher */ 2351 x[i] = 0.0; 2352 if (idx) { 2353 idx[i] = 0; /* in case ncols is zero */ 2354 for (j=0;j<ncols;j++) { /* find first implicit 0.0 in the row */ 2355 if (aj[j] > j) { 2356 idx[i] = j; 2357 break; 2358 } 2359 } 2360 } 2361 } 2362 for (j=0; j<ncols; j++){ 2363 if (PetscRealPart(x[i]) < PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2364 aa++; aj++; 2365 } 2366 } 2367 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2368 PetscFunctionReturn(0); 2369 } 2370 2371 #undef __FUNCT__ 2372 #define __FUNCT__ "MatGetRowMinAbs_SeqAIJ" 2373 PetscErrorCode MatGetRowMinAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2374 { 2375 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2376 PetscErrorCode ierr; 2377 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2378 PetscReal atmp; 2379 PetscScalar *x; 2380 MatScalar *aa; 2381 2382 PetscFunctionBegin; 2383 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2384 aa = a->a; 2385 ai = a->i; 2386 aj = a->j; 2387 2388 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2389 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2390 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2391 if (n != A->rmap->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2392 for (i=0; i<m; i++) { 2393 ncols = ai[1] - ai[0]; ai++; 2394 if (ncols) { 2395 /* Get first nonzero */ 2396 for(j = 0; j < ncols; j++) { 2397 atmp = PetscAbsScalar(aa[j]); 2398 if (atmp > 1.0e-12) {x[i] = atmp; if (idx) idx[i] = aj[j]; break;} 2399 } 2400 if (j == ncols) {x[i] = *aa; if (idx) idx[i] = *aj;} 2401 } else { 2402 x[i] = 0.0; if (idx) idx[i] = 0; 2403 } 2404 for(j = 0; j < ncols; j++) { 2405 atmp = PetscAbsScalar(*aa); 2406 if (atmp > 1.0e-12 && PetscAbsScalar(x[i]) > atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2407 aa++; aj++; 2408 } 2409 } 2410 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2411 PetscFunctionReturn(0); 2412 } 2413 2414 #undef __FUNCT__ 2415 #define __FUNCT__ "MatGetRowMin_SeqAIJ" 2416 PetscErrorCode MatGetRowMin_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2417 { 2418 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2419 PetscErrorCode ierr; 2420 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2421 PetscScalar *x; 2422 MatScalar *aa; 2423 2424 PetscFunctionBegin; 2425 if (A->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2426 aa = a->a; 2427 ai = a->i; 2428 aj = a->j; 2429 2430 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2431 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2432 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2433 if (n != A->rmap->n) SETERRQ(PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2434 for (i=0; i<m; i++) { 2435 ncols = ai[1] - ai[0]; ai++; 2436 if (ncols == A->cmap->n) { /* row is dense */ 2437 x[i] = *aa; if (idx) idx[i] = 0; 2438 } else { /* row is sparse so already KNOW minimum is 0.0 or lower */ 2439 x[i] = 0.0; 2440 if (idx) { /* find first implicit 0.0 in the row */ 2441 idx[i] = 0; /* in case ncols is zero */ 2442 for (j=0;j<ncols;j++) { 2443 if (aj[j] > j) { 2444 idx[i] = j; 2445 break; 2446 } 2447 } 2448 } 2449 } 2450 for (j=0; j<ncols; j++){ 2451 if (PetscRealPart(x[i]) > PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2452 aa++; aj++; 2453 } 2454 } 2455 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2456 PetscFunctionReturn(0); 2457 } 2458 extern PetscErrorCode PETSCMAT_DLLEXPORT MatFDColoringApply_AIJ(Mat,MatFDColoring,Vec,MatStructure*,void*); 2459 /* -------------------------------------------------------------------*/ 2460 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 2461 MatGetRow_SeqAIJ, 2462 MatRestoreRow_SeqAIJ, 2463 MatMult_SeqAIJ, 2464 /* 4*/ MatMultAdd_SeqAIJ, 2465 MatMultTranspose_SeqAIJ, 2466 MatMultTransposeAdd_SeqAIJ, 2467 0, 2468 0, 2469 0, 2470 /*10*/ 0, 2471 MatLUFactor_SeqAIJ, 2472 0, 2473 MatSOR_SeqAIJ, 2474 MatTranspose_SeqAIJ, 2475 /*15*/ MatGetInfo_SeqAIJ, 2476 MatEqual_SeqAIJ, 2477 MatGetDiagonal_SeqAIJ, 2478 MatDiagonalScale_SeqAIJ, 2479 MatNorm_SeqAIJ, 2480 /*20*/ 0, 2481 MatAssemblyEnd_SeqAIJ, 2482 MatSetOption_SeqAIJ, 2483 MatZeroEntries_SeqAIJ, 2484 /*24*/ MatZeroRows_SeqAIJ, 2485 0, 2486 0, 2487 0, 2488 0, 2489 /*29*/ MatSetUpPreallocation_SeqAIJ, 2490 0, 2491 0, 2492 MatGetArray_SeqAIJ, 2493 MatRestoreArray_SeqAIJ, 2494 /*34*/ MatDuplicate_SeqAIJ, 2495 0, 2496 0, 2497 MatILUFactor_SeqAIJ, 2498 0, 2499 /*39*/ MatAXPY_SeqAIJ, 2500 MatGetSubMatrices_SeqAIJ, 2501 MatIncreaseOverlap_SeqAIJ, 2502 MatGetValues_SeqAIJ, 2503 MatCopy_SeqAIJ, 2504 /*44*/ MatGetRowMax_SeqAIJ, 2505 MatScale_SeqAIJ, 2506 0, 2507 MatDiagonalSet_SeqAIJ, 2508 MatILUDTFactor_SeqAIJ, 2509 /*49*/ MatSetBlockSize_SeqAIJ, 2510 MatGetRowIJ_SeqAIJ, 2511 MatRestoreRowIJ_SeqAIJ, 2512 MatGetColumnIJ_SeqAIJ, 2513 MatRestoreColumnIJ_SeqAIJ, 2514 /*54*/ MatFDColoringCreate_SeqAIJ, 2515 0, 2516 0, 2517 MatPermute_SeqAIJ, 2518 0, 2519 /*59*/ 0, 2520 MatDestroy_SeqAIJ, 2521 MatView_SeqAIJ, 2522 0, 2523 0, 2524 /*64*/ 0, 2525 0, 2526 0, 2527 0, 2528 0, 2529 /*69*/ MatGetRowMaxAbs_SeqAIJ, 2530 MatGetRowMinAbs_SeqAIJ, 2531 0, 2532 MatSetColoring_SeqAIJ, 2533 #if defined(PETSC_HAVE_ADIC) 2534 MatSetValuesAdic_SeqAIJ, 2535 #else 2536 0, 2537 #endif 2538 /*74*/ MatSetValuesAdifor_SeqAIJ, 2539 MatFDColoringApply_AIJ, 2540 0, 2541 0, 2542 0, 2543 /*79*/ 0, 2544 0, 2545 0, 2546 0, 2547 MatLoad_SeqAIJ, 2548 /*84*/ MatIsSymmetric_SeqAIJ, 2549 MatIsHermitian_SeqAIJ, 2550 0, 2551 0, 2552 0, 2553 /*89*/ MatMatMult_SeqAIJ_SeqAIJ, 2554 MatMatMultSymbolic_SeqAIJ_SeqAIJ, 2555 MatMatMultNumeric_SeqAIJ_SeqAIJ, 2556 MatPtAP_Basic, 2557 MatPtAPSymbolic_SeqAIJ, 2558 /*94*/ MatPtAPNumeric_SeqAIJ, 2559 MatMatMultTranspose_SeqAIJ_SeqAIJ, 2560 MatMatMultTransposeSymbolic_SeqAIJ_SeqAIJ, 2561 MatMatMultTransposeNumeric_SeqAIJ_SeqAIJ, 2562 MatPtAPSymbolic_SeqAIJ_SeqAIJ, 2563 /*99*/ MatPtAPNumeric_SeqAIJ_SeqAIJ, 2564 0, 2565 0, 2566 MatConjugate_SeqAIJ, 2567 0, 2568 /*104*/MatSetValuesRow_SeqAIJ, 2569 MatRealPart_SeqAIJ, 2570 MatImaginaryPart_SeqAIJ, 2571 0, 2572 0, 2573 /*109*/0, 2574 0, 2575 MatGetRowMin_SeqAIJ, 2576 0, 2577 MatMissingDiagonal_SeqAIJ, 2578 /*114*/0, 2579 0, 2580 0, 2581 MatILUDTFactorSymbolic_SeqAIJ, 2582 MatILUDTFactorNumeric_SeqAIJ 2583 }; 2584 2585 EXTERN_C_BEGIN 2586 #undef __FUNCT__ 2587 #define __FUNCT__ "MatSeqAIJSetColumnIndices_SeqAIJ" 2588 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 2589 { 2590 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2591 PetscInt i,nz,n; 2592 2593 PetscFunctionBegin; 2594 2595 nz = aij->maxnz; 2596 n = mat->rmap->n; 2597 for (i=0; i<nz; i++) { 2598 aij->j[i] = indices[i]; 2599 } 2600 aij->nz = nz; 2601 for (i=0; i<n; i++) { 2602 aij->ilen[i] = aij->imax[i]; 2603 } 2604 2605 PetscFunctionReturn(0); 2606 } 2607 EXTERN_C_END 2608 2609 #undef __FUNCT__ 2610 #define __FUNCT__ "MatSeqAIJSetColumnIndices" 2611 /*@ 2612 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 2613 in the matrix. 2614 2615 Input Parameters: 2616 + mat - the SeqAIJ matrix 2617 - indices - the column indices 2618 2619 Level: advanced 2620 2621 Notes: 2622 This can be called if you have precomputed the nonzero structure of the 2623 matrix and want to provide it to the matrix object to improve the performance 2624 of the MatSetValues() operation. 2625 2626 You MUST have set the correct numbers of nonzeros per row in the call to 2627 MatCreateSeqAIJ(), and the columns indices MUST be sorted. 2628 2629 MUST be called before any calls to MatSetValues(); 2630 2631 The indices should start with zero, not one. 2632 2633 @*/ 2634 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 2635 { 2636 PetscErrorCode ierr,(*f)(Mat,PetscInt *); 2637 2638 PetscFunctionBegin; 2639 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 2640 PetscValidPointer(indices,2); 2641 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqAIJSetColumnIndices_C",(void (**)(void))&f);CHKERRQ(ierr); 2642 if (f) { 2643 ierr = (*f)(mat,indices);CHKERRQ(ierr); 2644 } else { 2645 SETERRQ(PETSC_ERR_SUP,"Wrong type of matrix to set column indices"); 2646 } 2647 PetscFunctionReturn(0); 2648 } 2649 2650 /* ----------------------------------------------------------------------------------------*/ 2651 2652 EXTERN_C_BEGIN 2653 #undef __FUNCT__ 2654 #define __FUNCT__ "MatStoreValues_SeqAIJ" 2655 PetscErrorCode PETSCMAT_DLLEXPORT MatStoreValues_SeqAIJ(Mat mat) 2656 { 2657 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2658 PetscErrorCode ierr; 2659 size_t nz = aij->i[mat->rmap->n]; 2660 2661 PetscFunctionBegin; 2662 if (aij->nonew != 1) { 2663 SETERRQ(PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 2664 } 2665 2666 /* allocate space for values if not already there */ 2667 if (!aij->saved_values) { 2668 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 2669 ierr = PetscLogObjectMemory(mat,(nz+1)*sizeof(PetscScalar));CHKERRQ(ierr); 2670 } 2671 2672 /* copy values over */ 2673 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2674 PetscFunctionReturn(0); 2675 } 2676 EXTERN_C_END 2677 2678 #undef __FUNCT__ 2679 #define __FUNCT__ "MatStoreValues" 2680 /*@ 2681 MatStoreValues - Stashes a copy of the matrix values; this allows, for 2682 example, reuse of the linear part of a Jacobian, while recomputing the 2683 nonlinear portion. 2684 2685 Collect on Mat 2686 2687 Input Parameters: 2688 . mat - the matrix (currently only AIJ matrices support this option) 2689 2690 Level: advanced 2691 2692 Common Usage, with SNESSolve(): 2693 $ Create Jacobian matrix 2694 $ Set linear terms into matrix 2695 $ Apply boundary conditions to matrix, at this time matrix must have 2696 $ final nonzero structure (i.e. setting the nonlinear terms and applying 2697 $ boundary conditions again will not change the nonzero structure 2698 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 2699 $ ierr = MatStoreValues(mat); 2700 $ Call SNESSetJacobian() with matrix 2701 $ In your Jacobian routine 2702 $ ierr = MatRetrieveValues(mat); 2703 $ Set nonlinear terms in matrix 2704 2705 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 2706 $ // build linear portion of Jacobian 2707 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 2708 $ ierr = MatStoreValues(mat); 2709 $ loop over nonlinear iterations 2710 $ ierr = MatRetrieveValues(mat); 2711 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 2712 $ // call MatAssemblyBegin/End() on matrix 2713 $ Solve linear system with Jacobian 2714 $ endloop 2715 2716 Notes: 2717 Matrix must already be assemblied before calling this routine 2718 Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before 2719 calling this routine. 2720 2721 When this is called multiple times it overwrites the previous set of stored values 2722 and does not allocated additional space. 2723 2724 .seealso: MatRetrieveValues() 2725 2726 @*/ 2727 PetscErrorCode PETSCMAT_DLLEXPORT MatStoreValues(Mat mat) 2728 { 2729 PetscErrorCode ierr,(*f)(Mat); 2730 2731 PetscFunctionBegin; 2732 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 2733 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2734 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2735 2736 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatStoreValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2737 if (f) { 2738 ierr = (*f)(mat);CHKERRQ(ierr); 2739 } else { 2740 SETERRQ(PETSC_ERR_SUP,"Wrong type of matrix to store values"); 2741 } 2742 PetscFunctionReturn(0); 2743 } 2744 2745 EXTERN_C_BEGIN 2746 #undef __FUNCT__ 2747 #define __FUNCT__ "MatRetrieveValues_SeqAIJ" 2748 PetscErrorCode PETSCMAT_DLLEXPORT MatRetrieveValues_SeqAIJ(Mat mat) 2749 { 2750 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2751 PetscErrorCode ierr; 2752 PetscInt nz = aij->i[mat->rmap->n]; 2753 2754 PetscFunctionBegin; 2755 if (aij->nonew != 1) { 2756 SETERRQ(PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 2757 } 2758 if (!aij->saved_values) { 2759 SETERRQ(PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 2760 } 2761 /* copy values over */ 2762 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 2763 PetscFunctionReturn(0); 2764 } 2765 EXTERN_C_END 2766 2767 #undef __FUNCT__ 2768 #define __FUNCT__ "MatRetrieveValues" 2769 /*@ 2770 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 2771 example, reuse of the linear part of a Jacobian, while recomputing the 2772 nonlinear portion. 2773 2774 Collect on Mat 2775 2776 Input Parameters: 2777 . mat - the matrix (currently on AIJ matrices support this option) 2778 2779 Level: advanced 2780 2781 .seealso: MatStoreValues() 2782 2783 @*/ 2784 PetscErrorCode PETSCMAT_DLLEXPORT MatRetrieveValues(Mat mat) 2785 { 2786 PetscErrorCode ierr,(*f)(Mat); 2787 2788 PetscFunctionBegin; 2789 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 2790 if (!mat->assembled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 2791 if (mat->factor) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2792 2793 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatRetrieveValues_C",(void (**)(void))&f);CHKERRQ(ierr); 2794 if (f) { 2795 ierr = (*f)(mat);CHKERRQ(ierr); 2796 } else { 2797 SETERRQ(PETSC_ERR_SUP,"Wrong type of matrix to retrieve values"); 2798 } 2799 PetscFunctionReturn(0); 2800 } 2801 2802 2803 /* --------------------------------------------------------------------------------*/ 2804 #undef __FUNCT__ 2805 #define __FUNCT__ "MatCreateSeqAIJ" 2806 /*@C 2807 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 2808 (the default parallel PETSc format). For good matrix assembly performance 2809 the user should preallocate the matrix storage by setting the parameter nz 2810 (or the array nnz). By setting these parameters accurately, performance 2811 during matrix assembly can be increased by more than a factor of 50. 2812 2813 Collective on MPI_Comm 2814 2815 Input Parameters: 2816 + comm - MPI communicator, set to PETSC_COMM_SELF 2817 . m - number of rows 2818 . n - number of columns 2819 . nz - number of nonzeros per row (same for all rows) 2820 - nnz - array containing the number of nonzeros in the various rows 2821 (possibly different for each row) or PETSC_NULL 2822 2823 Output Parameter: 2824 . A - the matrix 2825 2826 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 2827 MatXXXXSetPreallocation() paradgm instead of this routine directly. 2828 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 2829 2830 Notes: 2831 If nnz is given then nz is ignored 2832 2833 The AIJ format (also called the Yale sparse matrix format or 2834 compressed row storage), is fully compatible with standard Fortran 77 2835 storage. That is, the stored row and column indices can begin at 2836 either one (as in Fortran) or zero. See the users' manual for details. 2837 2838 Specify the preallocated storage with either nz or nnz (not both). 2839 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2840 allocation. For large problems you MUST preallocate memory or you 2841 will get TERRIBLE performance, see the users' manual chapter on matrices. 2842 2843 By default, this format uses inodes (identical nodes) when possible, to 2844 improve numerical efficiency of matrix-vector products and solves. We 2845 search for consecutive rows with the same nonzero structure, thereby 2846 reusing matrix information to achieve increased efficiency. 2847 2848 Options Database Keys: 2849 + -mat_no_inode - Do not use inodes 2850 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 2851 - -mat_aij_oneindex - Internally use indexing starting at 1 2852 rather than 0. Note that when calling MatSetValues(), 2853 the user still MUST index entries starting at 0! 2854 2855 Level: intermediate 2856 2857 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 2858 2859 @*/ 2860 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 2861 { 2862 PetscErrorCode ierr; 2863 2864 PetscFunctionBegin; 2865 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2866 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2867 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 2868 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,(PetscInt*)nnz);CHKERRQ(ierr); 2869 PetscFunctionReturn(0); 2870 } 2871 2872 #undef __FUNCT__ 2873 #define __FUNCT__ "MatSeqAIJSetPreallocation" 2874 /*@C 2875 MatSeqAIJSetPreallocation - For good matrix assembly performance 2876 the user should preallocate the matrix storage by setting the parameter nz 2877 (or the array nnz). By setting these parameters accurately, performance 2878 during matrix assembly can be increased by more than a factor of 50. 2879 2880 Collective on MPI_Comm 2881 2882 Input Parameters: 2883 + B - The matrix-free 2884 . nz - number of nonzeros per row (same for all rows) 2885 - nnz - array containing the number of nonzeros in the various rows 2886 (possibly different for each row) or PETSC_NULL 2887 2888 Notes: 2889 If nnz is given then nz is ignored 2890 2891 The AIJ format (also called the Yale sparse matrix format or 2892 compressed row storage), is fully compatible with standard Fortran 77 2893 storage. That is, the stored row and column indices can begin at 2894 either one (as in Fortran) or zero. See the users' manual for details. 2895 2896 Specify the preallocated storage with either nz or nnz (not both). 2897 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 2898 allocation. For large problems you MUST preallocate memory or you 2899 will get TERRIBLE performance, see the users' manual chapter on matrices. 2900 2901 You can call MatGetInfo() to get information on how effective the preallocation was; 2902 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 2903 You can also run with the option -info and look for messages with the string 2904 malloc in them to see if additional memory allocation was needed. 2905 2906 Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix 2907 entries or columns indices 2908 2909 By default, this format uses inodes (identical nodes) when possible, to 2910 improve numerical efficiency of matrix-vector products and solves. We 2911 search for consecutive rows with the same nonzero structure, thereby 2912 reusing matrix information to achieve increased efficiency. 2913 2914 Options Database Keys: 2915 + -mat_no_inode - Do not use inodes 2916 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 2917 - -mat_aij_oneindex - Internally use indexing starting at 1 2918 rather than 0. Note that when calling MatSetValues(), 2919 the user still MUST index entries starting at 0! 2920 2921 Level: intermediate 2922 2923 .seealso: MatCreate(), MatCreateMPIAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatGetInfo() 2924 2925 @*/ 2926 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 2927 { 2928 PetscErrorCode ierr,(*f)(Mat,PetscInt,const PetscInt[]); 2929 2930 PetscFunctionBegin; 2931 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 2932 if (f) { 2933 ierr = (*f)(B,nz,nnz);CHKERRQ(ierr); 2934 } 2935 PetscFunctionReturn(0); 2936 } 2937 2938 EXTERN_C_BEGIN 2939 #undef __FUNCT__ 2940 #define __FUNCT__ "MatSeqAIJSetPreallocation_SeqAIJ" 2941 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,PetscInt *nnz) 2942 { 2943 Mat_SeqAIJ *b; 2944 PetscTruth skipallocation = PETSC_FALSE; 2945 PetscErrorCode ierr; 2946 PetscInt i; 2947 2948 PetscFunctionBegin; 2949 2950 if (nz == MAT_SKIP_ALLOCATION) { 2951 skipallocation = PETSC_TRUE; 2952 nz = 0; 2953 } 2954 2955 ierr = PetscLayoutSetBlockSize(B->rmap,1);CHKERRQ(ierr); 2956 ierr = PetscLayoutSetBlockSize(B->cmap,1);CHKERRQ(ierr); 2957 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 2958 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 2959 2960 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 2961 if (nz < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 2962 if (nnz) { 2963 for (i=0; i<B->rmap->n; i++) { 2964 if (nnz[i] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %d value %d",i,nnz[i]); 2965 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); 2966 } 2967 } 2968 2969 B->preallocated = PETSC_TRUE; 2970 b = (Mat_SeqAIJ*)B->data; 2971 2972 if (!skipallocation) { 2973 if (!b->imax) { 2974 ierr = PetscMalloc2(B->rmap->n,PetscInt,&b->imax,B->rmap->n,PetscInt,&b->ilen);CHKERRQ(ierr); 2975 ierr = PetscLogObjectMemory(B,2*B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 2976 } 2977 if (!nnz) { 2978 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 2979 else if (nz <= 0) nz = 1; 2980 for (i=0; i<B->rmap->n; i++) b->imax[i] = nz; 2981 nz = nz*B->rmap->n; 2982 } else { 2983 nz = 0; 2984 for (i=0; i<B->rmap->n; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 2985 } 2986 /* b->ilen will count nonzeros in each row so far. */ 2987 for (i=0; i<B->rmap->n; i++) { b->ilen[i] = 0; } 2988 2989 /* allocate the matrix space */ 2990 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 2991 ierr = PetscMalloc3(nz,PetscScalar,&b->a,nz,PetscInt,&b->j,B->rmap->n+1,PetscInt,&b->i);CHKERRQ(ierr); 2992 ierr = PetscLogObjectMemory(B,(B->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 2993 b->i[0] = 0; 2994 for (i=1; i<B->rmap->n+1; i++) { 2995 b->i[i] = b->i[i-1] + b->imax[i-1]; 2996 } 2997 b->singlemalloc = PETSC_TRUE; 2998 b->free_a = PETSC_TRUE; 2999 b->free_ij = PETSC_TRUE; 3000 } else { 3001 b->free_a = PETSC_FALSE; 3002 b->free_ij = PETSC_FALSE; 3003 } 3004 3005 b->nz = 0; 3006 b->maxnz = nz; 3007 B->info.nz_unneeded = (double)b->maxnz; 3008 PetscFunctionReturn(0); 3009 } 3010 EXTERN_C_END 3011 3012 #undef __FUNCT__ 3013 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR" 3014 /*@ 3015 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 3016 3017 Input Parameters: 3018 + B - the matrix 3019 . i - the indices into j for the start of each row (starts with zero) 3020 . j - the column indices for each row (starts with zero) these must be sorted for each row 3021 - v - optional values in the matrix 3022 3023 Level: developer 3024 3025 The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays() 3026 3027 .keywords: matrix, aij, compressed row, sparse, sequential 3028 3029 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), SeqAIJ 3030 @*/ 3031 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 3032 { 3033 PetscErrorCode (*f)(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]); 3034 PetscErrorCode ierr; 3035 3036 PetscFunctionBegin; 3037 PetscValidHeaderSpecific(B,MAT_COOKIE,1); 3038 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C",(void (**)(void))&f);CHKERRQ(ierr); 3039 if (f) { 3040 ierr = (*f)(B,i,j,v);CHKERRQ(ierr); 3041 } 3042 PetscFunctionReturn(0); 3043 } 3044 3045 EXTERN_C_BEGIN 3046 #undef __FUNCT__ 3047 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR_SeqAIJ" 3048 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[]) 3049 { 3050 PetscInt i; 3051 PetscInt m,n; 3052 PetscInt nz; 3053 PetscInt *nnz, nz_max = 0; 3054 PetscScalar *values; 3055 PetscErrorCode ierr; 3056 3057 PetscFunctionBegin; 3058 ierr = MatGetSize(B, &m, &n);CHKERRQ(ierr); 3059 3060 if (Ii[0]) { 3061 SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %D", Ii[0]); 3062 } 3063 ierr = PetscMalloc((m+1) * sizeof(PetscInt), &nnz);CHKERRQ(ierr); 3064 for(i = 0; i < m; i++) { 3065 nz = Ii[i+1]- Ii[i]; 3066 nz_max = PetscMax(nz_max, nz); 3067 if (nz < 0) { 3068 SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE, "Local row %D has a negative number of columns %D", i, nnz); 3069 } 3070 nnz[i] = nz; 3071 } 3072 ierr = MatSeqAIJSetPreallocation(B, 0, nnz);CHKERRQ(ierr); 3073 ierr = PetscFree(nnz);CHKERRQ(ierr); 3074 3075 if (v) { 3076 values = (PetscScalar*) v; 3077 } else { 3078 ierr = PetscMalloc((nz_max+1)*sizeof(PetscScalar), &values);CHKERRQ(ierr); 3079 ierr = PetscMemzero(values, nz_max*sizeof(PetscScalar));CHKERRQ(ierr); 3080 } 3081 3082 for(i = 0; i < m; i++) { 3083 nz = Ii[i+1] - Ii[i]; 3084 ierr = MatSetValues_SeqAIJ(B, 1, &i, nz, J+Ii[i], values + (v ? Ii[i] : 0), INSERT_VALUES);CHKERRQ(ierr); 3085 } 3086 3087 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3088 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3089 3090 if (!v) { 3091 ierr = PetscFree(values);CHKERRQ(ierr); 3092 } 3093 PetscFunctionReturn(0); 3094 } 3095 EXTERN_C_END 3096 3097 #include "../src/mat/impls/dense/seq/dense.h" 3098 #include "private/petscaxpy.h" 3099 3100 #undef __FUNCT__ 3101 #define __FUNCT__ "MatMatMultNumeric_SeqDense_SeqAIJ" 3102 /* 3103 Computes (B'*A')' since computing B*A directly is untenable 3104 3105 n p p 3106 ( ) ( ) ( ) 3107 m ( A ) * n ( B ) = m ( C ) 3108 ( ) ( ) ( ) 3109 3110 */ 3111 PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C) 3112 { 3113 PetscErrorCode ierr; 3114 Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data; 3115 Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data; 3116 Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data; 3117 PetscInt i,n,m,q,p; 3118 const PetscInt *ii,*idx; 3119 const PetscScalar *b,*a,*a_q; 3120 PetscScalar *c,*c_q; 3121 3122 PetscFunctionBegin; 3123 m = A->rmap->n; 3124 n = A->cmap->n; 3125 p = B->cmap->n; 3126 a = sub_a->v; 3127 b = sub_b->a; 3128 c = sub_c->v; 3129 ierr = PetscMemzero(c,m*p*sizeof(PetscScalar));CHKERRQ(ierr); 3130 3131 ii = sub_b->i; 3132 idx = sub_b->j; 3133 for (i=0; i<n; i++) { 3134 q = ii[i+1] - ii[i]; 3135 while (q-->0) { 3136 c_q = c + m*(*idx); 3137 a_q = a + m*i; 3138 PetscAXPY(c_q,*b,a_q,m); 3139 idx++; 3140 b++; 3141 } 3142 } 3143 PetscFunctionReturn(0); 3144 } 3145 3146 #undef __FUNCT__ 3147 #define __FUNCT__ "MatMatMultSymbolic_SeqDense_SeqAIJ" 3148 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 3149 { 3150 PetscErrorCode ierr; 3151 PetscInt m=A->rmap->n,n=B->cmap->n; 3152 Mat Cmat; 3153 3154 PetscFunctionBegin; 3155 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); 3156 ierr = MatCreate(((PetscObject)A)->comm,&Cmat);CHKERRQ(ierr); 3157 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 3158 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 3159 ierr = MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);CHKERRQ(ierr); 3160 Cmat->assembled = PETSC_TRUE; 3161 *C = Cmat; 3162 PetscFunctionReturn(0); 3163 } 3164 3165 /* ----------------------------------------------------------------*/ 3166 #undef __FUNCT__ 3167 #define __FUNCT__ "MatMatMult_SeqDense_SeqAIJ" 3168 PetscErrorCode MatMatMult_SeqDense_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 3169 { 3170 PetscErrorCode ierr; 3171 3172 PetscFunctionBegin; 3173 if (scall == MAT_INITIAL_MATRIX){ 3174 ierr = MatMatMultSymbolic_SeqDense_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 3175 } 3176 ierr = MatMatMultNumeric_SeqDense_SeqAIJ(A,B,*C);CHKERRQ(ierr); 3177 PetscFunctionReturn(0); 3178 } 3179 3180 3181 /*MC 3182 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 3183 based on compressed sparse row format. 3184 3185 Options Database Keys: 3186 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 3187 3188 Level: beginner 3189 3190 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 3191 M*/ 3192 3193 EXTERN_C_BEGIN 3194 #if defined(PETSC_HAVE_PASTIX) 3195 extern PetscErrorCode MatGetFactor_seqaij_pastix(Mat,MatFactorType,Mat*); 3196 #endif 3197 #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SCALAR_SINGLE) && !defined(PETSC_USE_SCALAR_MAT_SINGLE) 3198 extern PetscErrorCode MatGetFactor_seqaij_essl(Mat,MatFactorType,Mat *); 3199 #endif 3200 extern PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqAIJ_SeqCRL(Mat,MatType,MatReuse,Mat*); 3201 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_petsc(Mat,MatFactorType,Mat*); 3202 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactorAvailable_seqaij_petsc(Mat,MatFactorType,PetscTruth *); 3203 #if defined(PETSC_HAVE_MUMPS) 3204 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_mumps(Mat,MatFactorType,Mat*); 3205 #endif 3206 #if defined(PETSC_HAVE_SUPERLU) 3207 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_superlu(Mat,MatFactorType,Mat*); 3208 #endif 3209 #if defined(PETSC_HAVE_SUPERLU_DIST) 3210 extern PetscErrorCode MatGetFactor_seqaij_superlu_dist(Mat,MatFactorType,Mat*); 3211 #endif 3212 #if defined(PETSC_HAVE_SPOOLES) 3213 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_spooles(Mat,MatFactorType,Mat*); 3214 #endif 3215 #if defined(PETSC_HAVE_UMFPACK) 3216 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_umfpack(Mat,MatFactorType,Mat*); 3217 #endif 3218 #if defined(PETSC_HAVE_LUSOL) 3219 extern PetscErrorCode PETSCMAT_DLLEXPORT MatGetFactor_seqaij_lusol(Mat,MatFactorType,Mat*); 3220 #endif 3221 EXTERN_C_END 3222 3223 3224 EXTERN_C_BEGIN 3225 #undef __FUNCT__ 3226 #define __FUNCT__ "MatCreate_SeqAIJ" 3227 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_SeqAIJ(Mat B) 3228 { 3229 Mat_SeqAIJ *b; 3230 PetscErrorCode ierr; 3231 PetscMPIInt size; 3232 3233 PetscFunctionBegin; 3234 ierr = MPI_Comm_size(((PetscObject)B)->comm,&size);CHKERRQ(ierr); 3235 if (size > 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 3236 3237 ierr = PetscNewLog(B,Mat_SeqAIJ,&b);CHKERRQ(ierr); 3238 B->data = (void*)b; 3239 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 3240 B->mapping = 0; 3241 b->row = 0; 3242 b->col = 0; 3243 b->icol = 0; 3244 b->reallocs = 0; 3245 b->ignorezeroentries = PETSC_FALSE; 3246 b->roworiented = PETSC_TRUE; 3247 b->nonew = 0; 3248 b->diag = 0; 3249 b->solve_work = 0; 3250 B->spptr = 0; 3251 b->saved_values = 0; 3252 b->idiag = 0; 3253 b->mdiag = 0; 3254 b->ssor_work = 0; 3255 b->omega = 1.0; 3256 b->fshift = 0.0; 3257 b->idiagvalid = PETSC_FALSE; 3258 b->keepnonzeropattern = PETSC_FALSE; 3259 b->xtoy = 0; 3260 b->XtoY = 0; 3261 b->compressedrow.use = PETSC_FALSE; 3262 b->compressedrow.nrows = B->rmap->n; 3263 b->compressedrow.i = PETSC_NULL; 3264 b->compressedrow.rindex = PETSC_NULL; 3265 b->compressedrow.checked = PETSC_FALSE; 3266 B->same_nonzero = PETSC_FALSE; 3267 3268 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3269 #if defined(PETSC_HAVE_PASTIX) 3270 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_pastix_C", 3271 "MatGetFactor_seqaij_pastix", 3272 MatGetFactor_seqaij_pastix);CHKERRQ(ierr); 3273 #endif 3274 #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SCALAR_SINGLE) && !defined(PETSC_USE_SCALAR_MAT_SINGLE) 3275 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_essl_C", 3276 "MatGetFactor_seqaij_essl", 3277 MatGetFactor_seqaij_essl);CHKERRQ(ierr); 3278 #endif 3279 #if defined(PETSC_HAVE_SUPERLU) 3280 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_superlu_C", 3281 "MatGetFactor_seqaij_superlu", 3282 MatGetFactor_seqaij_superlu);CHKERRQ(ierr); 3283 #endif 3284 #if defined(PETSC_HAVE_SUPERLU_DIST) 3285 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_superlu_dist_C", 3286 "MatGetFactor_seqaij_superlu_dist", 3287 MatGetFactor_seqaij_superlu_dist);CHKERRQ(ierr); 3288 #endif 3289 #if defined(PETSC_HAVE_SPOOLES) 3290 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_spooles_C", 3291 "MatGetFactor_seqaij_spooles", 3292 MatGetFactor_seqaij_spooles);CHKERRQ(ierr); 3293 #endif 3294 #if defined(PETSC_HAVE_MUMPS) 3295 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_mumps_C", 3296 "MatGetFactor_seqaij_mumps", 3297 MatGetFactor_seqaij_mumps);CHKERRQ(ierr); 3298 #endif 3299 #if defined(PETSC_HAVE_UMFPACK) 3300 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_umfpack_C", 3301 "MatGetFactor_seqaij_umfpack", 3302 MatGetFactor_seqaij_umfpack);CHKERRQ(ierr); 3303 #endif 3304 #if defined(PETSC_HAVE_LUSOL) 3305 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_lusol_C", 3306 "MatGetFactor_seqaij_lusol", 3307 MatGetFactor_seqaij_lusol);CHKERRQ(ierr); 3308 #endif 3309 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_petsc_C", 3310 "MatGetFactor_seqaij_petsc", 3311 MatGetFactor_seqaij_petsc);CHKERRQ(ierr); 3312 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactorAvailable_petsc_C", 3313 "MatGetFactorAvailable_seqaij_petsc", 3314 MatGetFactorAvailable_seqaij_petsc);CHKERRQ(ierr); 3315 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetColumnIndices_C", 3316 "MatSeqAIJSetColumnIndices_SeqAIJ", 3317 MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 3318 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 3319 "MatStoreValues_SeqAIJ", 3320 MatStoreValues_SeqAIJ);CHKERRQ(ierr); 3321 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 3322 "MatRetrieveValues_SeqAIJ", 3323 MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 3324 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqsbaij_C", 3325 "MatConvert_SeqAIJ_SeqSBAIJ", 3326 MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 3327 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqbaij_C", 3328 "MatConvert_SeqAIJ_SeqBAIJ", 3329 MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 3330 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqcsrperm_C", 3331 "MatConvert_SeqAIJ_SeqCSRPERM", 3332 MatConvert_SeqAIJ_SeqCSRPERM);CHKERRQ(ierr); 3333 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqcrl_C", 3334 "MatConvert_SeqAIJ_SeqCRL", 3335 MatConvert_SeqAIJ_SeqCRL);CHKERRQ(ierr); 3336 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsTranspose_C", 3337 "MatIsTranspose_SeqAIJ", 3338 MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3339 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsHermitianTranspose_C", 3340 "MatIsHermitianTranspose_SeqAIJ", 3341 MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3342 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocation_C", 3343 "MatSeqAIJSetPreallocation_SeqAIJ", 3344 MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 3345 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C", 3346 "MatSeqAIJSetPreallocationCSR_SeqAIJ", 3347 MatSeqAIJSetPreallocationCSR_SeqAIJ);CHKERRQ(ierr); 3348 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatReorderForNonzeroDiagonal_C", 3349 "MatReorderForNonzeroDiagonal_SeqAIJ", 3350 MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 3351 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMult_seqdense_seqaij_C", 3352 "MatMatMult_SeqDense_SeqAIJ", 3353 MatMatMult_SeqDense_SeqAIJ);CHKERRQ(ierr); 3354 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultSymbolic_seqdense_seqaij_C", 3355 "MatMatMultSymbolic_SeqDense_SeqAIJ", 3356 MatMatMultSymbolic_SeqDense_SeqAIJ);CHKERRQ(ierr); 3357 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultNumeric_seqdense_seqaij_C", 3358 "MatMatMultNumeric_SeqDense_SeqAIJ", 3359 MatMatMultNumeric_SeqDense_SeqAIJ);CHKERRQ(ierr); 3360 ierr = MatCreate_SeqAIJ_Inode(B);CHKERRQ(ierr); 3361 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3362 PetscFunctionReturn(0); 3363 } 3364 EXTERN_C_END 3365 3366 #undef __FUNCT__ 3367 #define __FUNCT__ "MatDuplicateNoCreate_SeqAIJ" 3368 /* 3369 Given a matrix generated with MatGetFactor() duplicates all the information in A into B 3370 */ 3371 PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscTruth mallocmatspace) 3372 { 3373 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 3374 PetscErrorCode ierr; 3375 PetscInt i,m = A->rmap->n; 3376 3377 PetscFunctionBegin; 3378 c = (Mat_SeqAIJ*)C->data; 3379 3380 C->factor = A->factor; 3381 3382 c->row = 0; 3383 c->col = 0; 3384 c->icol = 0; 3385 c->reallocs = 0; 3386 3387 C->assembled = PETSC_TRUE; 3388 3389 ierr = PetscLayoutSetBlockSize(C->rmap,1);CHKERRQ(ierr); 3390 ierr = PetscLayoutSetBlockSize(C->cmap,1);CHKERRQ(ierr); 3391 ierr = PetscLayoutSetUp(C->rmap);CHKERRQ(ierr); 3392 ierr = PetscLayoutSetUp(C->cmap);CHKERRQ(ierr); 3393 3394 ierr = PetscMalloc2(m,PetscInt,&c->imax,m,PetscInt,&c->ilen);CHKERRQ(ierr); 3395 ierr = PetscLogObjectMemory(C, 2*m*sizeof(PetscInt));CHKERRQ(ierr); 3396 for (i=0; i<m; i++) { 3397 c->imax[i] = a->imax[i]; 3398 c->ilen[i] = a->ilen[i]; 3399 } 3400 3401 /* allocate the matrix space */ 3402 if (mallocmatspace){ 3403 ierr = PetscMalloc3(a->i[m],PetscScalar,&c->a,a->i[m],PetscInt,&c->j,m+1,PetscInt,&c->i);CHKERRQ(ierr); 3404 ierr = PetscLogObjectMemory(C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3405 c->singlemalloc = PETSC_TRUE; 3406 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3407 if (m > 0) { 3408 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(PetscInt));CHKERRQ(ierr); 3409 if (cpvalues == MAT_COPY_VALUES) { 3410 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 3411 } else { 3412 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 3413 } 3414 } 3415 } 3416 3417 c->ignorezeroentries = a->ignorezeroentries; 3418 c->roworiented = a->roworiented; 3419 c->nonew = a->nonew; 3420 if (a->diag) { 3421 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 3422 ierr = PetscLogObjectMemory(C,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3423 for (i=0; i<m; i++) { 3424 c->diag[i] = a->diag[i]; 3425 } 3426 } else c->diag = 0; 3427 c->solve_work = 0; 3428 c->saved_values = 0; 3429 c->idiag = 0; 3430 c->ssor_work = 0; 3431 c->keepnonzeropattern = a->keepnonzeropattern; 3432 c->free_a = PETSC_TRUE; 3433 c->free_ij = PETSC_TRUE; 3434 c->xtoy = 0; 3435 c->XtoY = 0; 3436 3437 c->nz = a->nz; 3438 c->maxnz = a->maxnz; 3439 C->preallocated = PETSC_TRUE; 3440 3441 c->compressedrow.use = a->compressedrow.use; 3442 c->compressedrow.nrows = a->compressedrow.nrows; 3443 c->compressedrow.checked = a->compressedrow.checked; 3444 if ( a->compressedrow.checked && a->compressedrow.use){ 3445 i = a->compressedrow.nrows; 3446 ierr = PetscMalloc((2*i+1)*sizeof(PetscInt),&c->compressedrow.i);CHKERRQ(ierr); 3447 c->compressedrow.rindex = c->compressedrow.i + i + 1; 3448 ierr = PetscMemcpy(c->compressedrow.i,a->compressedrow.i,(i+1)*sizeof(PetscInt));CHKERRQ(ierr); 3449 ierr = PetscMemcpy(c->compressedrow.rindex,a->compressedrow.rindex,i*sizeof(PetscInt));CHKERRQ(ierr); 3450 } else { 3451 c->compressedrow.use = PETSC_FALSE; 3452 c->compressedrow.i = PETSC_NULL; 3453 c->compressedrow.rindex = PETSC_NULL; 3454 } 3455 C->same_nonzero = A->same_nonzero; 3456 ierr = MatDuplicate_SeqAIJ_Inode(A,cpvalues,&C);CHKERRQ(ierr); 3457 3458 ierr = PetscFListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 3459 PetscFunctionReturn(0); 3460 } 3461 3462 #undef __FUNCT__ 3463 #define __FUNCT__ "MatDuplicate_SeqAIJ" 3464 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 3465 { 3466 PetscErrorCode ierr; 3467 PetscInt n = A->rmap->n; 3468 3469 PetscFunctionBegin; 3470 ierr = MatCreate(((PetscObject)A)->comm,B);CHKERRQ(ierr); 3471 ierr = MatSetSizes(*B,n,n,n,n);CHKERRQ(ierr); 3472 ierr = MatSetType(*B,MATSEQAIJ);CHKERRQ(ierr); 3473 ierr = MatDuplicateNoCreate_SeqAIJ(*B,A,cpvalues,PETSC_TRUE);CHKERRQ(ierr); 3474 PetscFunctionReturn(0); 3475 } 3476 3477 #undef __FUNCT__ 3478 #define __FUNCT__ "MatLoad_SeqAIJ" 3479 PetscErrorCode MatLoad_SeqAIJ(PetscViewer viewer, const MatType type,Mat *A) 3480 { 3481 Mat_SeqAIJ *a; 3482 Mat B; 3483 PetscErrorCode ierr; 3484 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N; 3485 int fd; 3486 PetscMPIInt size; 3487 MPI_Comm comm; 3488 3489 PetscFunctionBegin; 3490 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 3491 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 3492 if (size > 1) SETERRQ(PETSC_ERR_ARG_SIZ,"view must have one processor"); 3493 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 3494 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 3495 if (header[0] != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 3496 M = header[1]; N = header[2]; nz = header[3]; 3497 3498 if (nz < 0) { 3499 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 3500 } 3501 3502 /* read in row lengths */ 3503 ierr = PetscMalloc(M*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 3504 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 3505 3506 /* check if sum of rowlengths is same as nz */ 3507 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 3508 if (sum != nz) SETERRQ2(PETSC_ERR_FILE_READ,"Inconsistant matrix data in file. no-nonzeros = %d, sum-row-lengths = %d\n",nz,sum); 3509 3510 /* create our matrix */ 3511 ierr = MatCreate(comm,&B);CHKERRQ(ierr); 3512 ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 3513 ierr = MatSetType(B,type);CHKERRQ(ierr); 3514 ierr = MatSeqAIJSetPreallocation_SeqAIJ(B,0,rowlengths);CHKERRQ(ierr); 3515 a = (Mat_SeqAIJ*)B->data; 3516 3517 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 3518 3519 /* read in nonzero values */ 3520 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 3521 3522 /* set matrix "i" values */ 3523 a->i[0] = 0; 3524 for (i=1; i<= M; i++) { 3525 a->i[i] = a->i[i-1] + rowlengths[i-1]; 3526 a->ilen[i-1] = rowlengths[i-1]; 3527 } 3528 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 3529 3530 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3531 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3532 *A = B; 3533 PetscFunctionReturn(0); 3534 } 3535 3536 #undef __FUNCT__ 3537 #define __FUNCT__ "MatEqual_SeqAIJ" 3538 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscTruth* flg) 3539 { 3540 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data,*b = (Mat_SeqAIJ *)B->data; 3541 PetscErrorCode ierr; 3542 #if defined(PETSC_USE_COMPLEX) 3543 PetscInt k; 3544 #endif 3545 3546 PetscFunctionBegin; 3547 /* If the matrix dimensions are not equal,or no of nonzeros */ 3548 if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) ||(a->nz != b->nz)) { 3549 *flg = PETSC_FALSE; 3550 PetscFunctionReturn(0); 3551 } 3552 3553 /* if the a->i are the same */ 3554 ierr = PetscMemcmp(a->i,b->i,(A->rmap->n+1)*sizeof(PetscInt),flg);CHKERRQ(ierr); 3555 if (!*flg) PetscFunctionReturn(0); 3556 3557 /* if a->j are the same */ 3558 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),flg);CHKERRQ(ierr); 3559 if (!*flg) PetscFunctionReturn(0); 3560 3561 /* if a->a are the same */ 3562 #if defined(PETSC_USE_COMPLEX) 3563 for (k=0; k<a->nz; k++){ 3564 if (PetscRealPart(a->a[k]) != PetscRealPart(b->a[k]) || PetscImaginaryPart(a->a[k]) != PetscImaginaryPart(b->a[k])){ 3565 *flg = PETSC_FALSE; 3566 PetscFunctionReturn(0); 3567 } 3568 } 3569 #else 3570 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 3571 #endif 3572 PetscFunctionReturn(0); 3573 } 3574 3575 #undef __FUNCT__ 3576 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 3577 /*@ 3578 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 3579 provided by the user. 3580 3581 Collective on MPI_Comm 3582 3583 Input Parameters: 3584 + comm - must be an MPI communicator of size 1 3585 . m - number of rows 3586 . n - number of columns 3587 . i - row indices 3588 . j - column indices 3589 - a - matrix values 3590 3591 Output Parameter: 3592 . mat - the matrix 3593 3594 Level: intermediate 3595 3596 Notes: 3597 The i, j, and a arrays are not copied by this routine, the user must free these arrays 3598 once the matrix is destroyed 3599 3600 You cannot set new nonzero locations into this matrix, that will generate an error. 3601 3602 The i and j indices are 0 based 3603 3604 The format which is used for the sparse matrix input, is equivalent to a 3605 row-major ordering.. i.e for the following matrix, the input data expected is 3606 as shown: 3607 3608 1 0 0 3609 2 0 3 3610 4 5 6 3611 3612 i = {0,1,3,6} [size = nrow+1 = 3+1] 3613 j = {0,0,2,0,1,2} [size = nz = 6]; values must be sorted for each row 3614 v = {1,2,3,4,5,6} [size = nz = 6] 3615 3616 3617 .seealso: MatCreate(), MatCreateMPIAIJ(), MatCreateSeqAIJ(), MatCreateMPIAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 3618 3619 @*/ 3620 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt* i,PetscInt*j,PetscScalar *a,Mat *mat) 3621 { 3622 PetscErrorCode ierr; 3623 PetscInt ii; 3624 Mat_SeqAIJ *aij; 3625 #if defined(PETSC_USE_DEBUG) 3626 PetscInt jj; 3627 #endif 3628 3629 PetscFunctionBegin; 3630 if (i[0]) { 3631 SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 3632 } 3633 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 3634 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 3635 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 3636 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 3637 aij = (Mat_SeqAIJ*)(*mat)->data; 3638 ierr = PetscMalloc2(m,PetscInt,&aij->imax,m,PetscInt,&aij->ilen);CHKERRQ(ierr); 3639 3640 aij->i = i; 3641 aij->j = j; 3642 aij->a = a; 3643 aij->singlemalloc = PETSC_FALSE; 3644 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 3645 aij->free_a = PETSC_FALSE; 3646 aij->free_ij = PETSC_FALSE; 3647 3648 for (ii=0; ii<m; ii++) { 3649 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 3650 #if defined(PETSC_USE_DEBUG) 3651 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]); 3652 for (jj=i[ii]+1; jj<i[ii+1]; jj++) { 3653 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); 3654 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); 3655 } 3656 #endif 3657 } 3658 #if defined(PETSC_USE_DEBUG) 3659 for (ii=0; ii<aij->i[m]; ii++) { 3660 if (j[ii] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 3661 if (j[ii] > n - 1) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column index to large at location = %d index = %d",ii,j[ii]); 3662 } 3663 #endif 3664 3665 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3666 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3667 PetscFunctionReturn(0); 3668 } 3669 3670 #undef __FUNCT__ 3671 #define __FUNCT__ "MatSetColoring_SeqAIJ" 3672 PetscErrorCode MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 3673 { 3674 PetscErrorCode ierr; 3675 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3676 3677 PetscFunctionBegin; 3678 if (coloring->ctype == IS_COLORING_GLOBAL) { 3679 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 3680 a->coloring = coloring; 3681 } else if (coloring->ctype == IS_COLORING_GHOSTED) { 3682 PetscInt i,*larray; 3683 ISColoring ocoloring; 3684 ISColoringValue *colors; 3685 3686 /* set coloring for diagonal portion */ 3687 ierr = PetscMalloc((A->cmap->n+1)*sizeof(PetscInt),&larray);CHKERRQ(ierr); 3688 for (i=0; i<A->cmap->n; i++) { 3689 larray[i] = i; 3690 } 3691 ierr = ISGlobalToLocalMappingApply(A->mapping,IS_GTOLM_MASK,A->cmap->n,larray,PETSC_NULL,larray);CHKERRQ(ierr); 3692 ierr = PetscMalloc((A->cmap->n+1)*sizeof(ISColoringValue),&colors);CHKERRQ(ierr); 3693 for (i=0; i<A->cmap->n; i++) { 3694 colors[i] = coloring->colors[larray[i]]; 3695 } 3696 ierr = PetscFree(larray);CHKERRQ(ierr); 3697 ierr = ISColoringCreate(PETSC_COMM_SELF,coloring->n,A->cmap->n,colors,&ocoloring);CHKERRQ(ierr); 3698 a->coloring = ocoloring; 3699 } 3700 PetscFunctionReturn(0); 3701 } 3702 3703 #if defined(PETSC_HAVE_ADIC) 3704 EXTERN_C_BEGIN 3705 #include "adic/ad_utils.h" 3706 EXTERN_C_END 3707 3708 #undef __FUNCT__ 3709 #define __FUNCT__ "MatSetValuesAdic_SeqAIJ" 3710 PetscErrorCode MatSetValuesAdic_SeqAIJ(Mat A,void *advalues) 3711 { 3712 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3713 PetscInt m = A->rmap->n,*ii = a->i,*jj = a->j,nz,i,j,nlen; 3714 PetscScalar *v = a->a,*values = ((PetscScalar*)advalues)+1; 3715 ISColoringValue *color; 3716 3717 PetscFunctionBegin; 3718 if (!a->coloring) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 3719 nlen = PetscADGetDerivTypeSize()/sizeof(PetscScalar); 3720 color = a->coloring->colors; 3721 /* loop over rows */ 3722 for (i=0; i<m; i++) { 3723 nz = ii[i+1] - ii[i]; 3724 /* loop over columns putting computed value into matrix */ 3725 for (j=0; j<nz; j++) { 3726 *v++ = values[color[*jj++]]; 3727 } 3728 values += nlen; /* jump to next row of derivatives */ 3729 } 3730 PetscFunctionReturn(0); 3731 } 3732 #endif 3733 3734 #undef __FUNCT__ 3735 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 3736 PetscErrorCode MatSetValuesAdifor_SeqAIJ(Mat A,PetscInt nl,void *advalues) 3737 { 3738 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3739 PetscInt m = A->rmap->n,*ii = a->i,*jj = a->j,nz,i,j; 3740 MatScalar *v = a->a; 3741 PetscScalar *values = (PetscScalar *)advalues; 3742 ISColoringValue *color; 3743 3744 PetscFunctionBegin; 3745 if (!a->coloring) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 3746 color = a->coloring->colors; 3747 /* loop over rows */ 3748 for (i=0; i<m; i++) { 3749 nz = ii[i+1] - ii[i]; 3750 /* loop over columns putting computed value into matrix */ 3751 for (j=0; j<nz; j++) { 3752 *v++ = values[color[*jj++]]; 3753 } 3754 values += nl; /* jump to next row of derivatives */ 3755 } 3756 PetscFunctionReturn(0); 3757 } 3758 3759 /* 3760 Special version for direct calls from Fortran 3761 */ 3762 #if defined(PETSC_HAVE_FORTRAN_CAPS) 3763 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 3764 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 3765 #define matsetvaluesseqaij_ matsetvaluesseqaij 3766 #endif 3767 3768 /* Change these macros so can be used in void function */ 3769 #undef CHKERRQ 3770 #define CHKERRQ(ierr) CHKERRABORT(((PetscObject)A)->comm,ierr) 3771 #undef SETERRQ2 3772 #define SETERRQ2(ierr,b,c,d) CHKERRABORT(((PetscObject)A)->comm,ierr) 3773 3774 EXTERN_C_BEGIN 3775 #undef __FUNCT__ 3776 #define __FUNCT__ "matsetvaluesseqaij_" 3777 void PETSC_STDCALL matsetvaluesseqaij_(Mat *AA,PetscInt *mm,const PetscInt im[],PetscInt *nn,const PetscInt in[],const PetscScalar v[],InsertMode *isis, PetscErrorCode *_ierr) 3778 { 3779 Mat A = *AA; 3780 PetscInt m = *mm, n = *nn; 3781 InsertMode is = *isis; 3782 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3783 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 3784 PetscInt *imax,*ai,*ailen; 3785 PetscErrorCode ierr; 3786 PetscInt *aj,nonew = a->nonew,lastcol = -1; 3787 MatScalar *ap,value,*aa; 3788 PetscTruth ignorezeroentries = a->ignorezeroentries; 3789 PetscTruth roworiented = a->roworiented; 3790 3791 PetscFunctionBegin; 3792 ierr = MatPreallocated(A);CHKERRQ(ierr); 3793 imax = a->imax; 3794 ai = a->i; 3795 ailen = a->ilen; 3796 aj = a->j; 3797 aa = a->a; 3798 3799 for (k=0; k<m; k++) { /* loop over added rows */ 3800 row = im[k]; 3801 if (row < 0) continue; 3802 #if defined(PETSC_USE_DEBUG) 3803 if (row >= A->rmap->n) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 3804 #endif 3805 rp = aj + ai[row]; ap = aa + ai[row]; 3806 rmax = imax[row]; nrow = ailen[row]; 3807 low = 0; 3808 high = nrow; 3809 for (l=0; l<n; l++) { /* loop over added columns */ 3810 if (in[l] < 0) continue; 3811 #if defined(PETSC_USE_DEBUG) 3812 if (in[l] >= A->cmap->n) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 3813 #endif 3814 col = in[l]; 3815 if (roworiented) { 3816 value = v[l + k*n]; 3817 } else { 3818 value = v[k + l*m]; 3819 } 3820 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 3821 3822 if (col <= lastcol) low = 0; else high = nrow; 3823 lastcol = col; 3824 while (high-low > 5) { 3825 t = (low+high)/2; 3826 if (rp[t] > col) high = t; 3827 else low = t; 3828 } 3829 for (i=low; i<high; i++) { 3830 if (rp[i] > col) break; 3831 if (rp[i] == col) { 3832 if (is == ADD_VALUES) ap[i] += value; 3833 else ap[i] = value; 3834 goto noinsert; 3835 } 3836 } 3837 if (value == 0.0 && ignorezeroentries) goto noinsert; 3838 if (nonew == 1) goto noinsert; 3839 if (nonew == -1) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 3840 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 3841 N = nrow++ - 1; a->nz++; high++; 3842 /* shift up all the later entries in this row */ 3843 for (ii=N; ii>=i; ii--) { 3844 rp[ii+1] = rp[ii]; 3845 ap[ii+1] = ap[ii]; 3846 } 3847 rp[i] = col; 3848 ap[i] = value; 3849 noinsert:; 3850 low = i + 1; 3851 } 3852 ailen[row] = nrow; 3853 } 3854 A->same_nonzero = PETSC_FALSE; 3855 PetscFunctionReturnVoid(); 3856 } 3857 EXTERN_C_END 3858