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