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