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