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