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