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