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