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