Actual source code: aijbaij.c

petsc-3.14.6 2021-03-30
Report Typos and Errors

  2: #include <../src/mat/impls/baij/seq/baij.h>

  4: PETSC_INTERN PetscErrorCode MatConvert_SeqBAIJ_SeqAIJ(Mat A, MatType newtype,MatReuse reuse,Mat *newmat)
  5: {
  6:   Mat            B;
  7:   Mat_SeqAIJ     *b;
  8:   PetscBool      roworiented;
  9:   Mat_SeqBAIJ    *a = (Mat_SeqBAIJ*)A->data;
 11:   PetscInt       bs = A->rmap->bs,*ai = a->i,*aj = a->j,n = A->rmap->N/bs,i,j,k;
 12:   PetscInt       *rowlengths,*rows,*cols,maxlen = 0,ncols;
 13:   MatScalar      *aa = a->a;

 16:   if (reuse == MAT_REUSE_MATRIX) {
 17:     B = *newmat;
 18:     for (i=0; i<n; i++) {
 19:       maxlen = PetscMax(maxlen,(ai[i+1] - ai[i]));
 20:     }
 21:   } else {
 22:     PetscMalloc1(n*bs,&rowlengths);
 23:     for (i=0; i<n; i++) {
 24:       maxlen = PetscMax(maxlen,(ai[i+1] - ai[i]));
 25:       for (j=0; j<bs; j++) {
 26:         rowlengths[i*bs+j] = bs*(ai[i+1] - ai[i]);
 27:       }
 28:     }
 29:     MatCreate(PetscObjectComm((PetscObject)A),&B);
 30:     MatSetType(B,MATSEQAIJ);
 31:     MatSetSizes(B,A->rmap->n,A->cmap->n,A->rmap->N,A->cmap->N);
 32:     MatSetBlockSizes(B,A->rmap->bs,A->cmap->bs);
 33:     MatSeqAIJSetPreallocation(B,0,rowlengths);
 34:     PetscFree(rowlengths);
 35:   }
 36:   b = (Mat_SeqAIJ*)B->data;
 37:   roworiented = b->roworiented;

 39:   MatSetOption(B,MAT_ROW_ORIENTED,PETSC_FALSE);
 40:   PetscMalloc1(bs,&rows);
 41:   PetscMalloc1(bs*maxlen,&cols);
 42:   for (i=0; i<n; i++) {
 43:     for (j=0; j<bs; j++) {
 44:       rows[j] = i*bs+j;
 45:     }
 46:     ncols = ai[i+1] - ai[i];
 47:     for (k=0; k<ncols; k++) {
 48:       for (j=0; j<bs; j++) {
 49:         cols[k*bs+j] = bs*(*aj) + j;
 50:       }
 51:       aj++;
 52:     }
 53:     MatSetValues(B,bs,rows,bs*ncols,cols,aa,INSERT_VALUES);
 54:     aa  += ncols*bs*bs;
 55:   }
 56:   PetscFree(cols);
 57:   PetscFree(rows);
 58:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
 59:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
 60:   MatSetOption(B,MAT_ROW_ORIENTED,roworiented);

 62:   if (reuse == MAT_INPLACE_MATRIX) {
 63:     MatHeaderReplace(A,&B);
 64:   } else *newmat = B;
 65:   return(0);
 66: }

 68: #include <../src/mat/impls/aij/seq/aij.h>

 70: PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqBAIJ(Mat A,MatType newtype,MatReuse reuse,Mat *newmat)
 71: {
 72:   Mat            B;
 73:   Mat_SeqAIJ     *a = (Mat_SeqAIJ*)A->data;
 74:   Mat_SeqBAIJ    *b;
 76:   PetscInt       *ai=a->i,m=A->rmap->N,n=A->cmap->N,i,*rowlengths,bs=PetscAbs(A->rmap->bs);

 79:   if (reuse != MAT_REUSE_MATRIX) {
 80:     PetscMalloc1(m/bs,&rowlengths);
 81:     for (i=0; i<m/bs; i++) {
 82:       rowlengths[i] = (ai[i*bs+1] - ai[i*bs])/bs;
 83:     }
 84:     MatCreate(PetscObjectComm((PetscObject)A),&B);
 85:     MatSetSizes(B,m,n,m,n);
 86:     MatSetType(B,MATSEQBAIJ);
 87:     MatSeqBAIJSetPreallocation(B,bs,0,rowlengths);
 88:     PetscFree(rowlengths);
 89:   } else B = *newmat;

 91:   if (bs == 1) {
 92:     b = (Mat_SeqBAIJ*)(B->data);

 94:     PetscArraycpy(b->i,a->i,m+1);
 95:     PetscArraycpy(b->ilen,a->ilen,m);
 96:     PetscArraycpy(b->j,a->j,a->nz);
 97:     PetscArraycpy(b->a,a->a,a->nz);

 99:     MatSetOption(B,MAT_ROW_ORIENTED,PETSC_TRUE);
100:     MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
101:     MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);
102:   } else {
103:     /* reuse may not be equal to MAT_REUSE_MATRIX, but the basic converter will reallocate or replace newmat if this value is not used */
104:     /* if reuse is equal to MAT_INITIAL_MATRIX, it has been appropriately preallocated before                                          */
105:     /*                      MAT_INPLACE_MATRIX, it will be replaced with MatHeaderReplace below                                        */
106:     MatConvert_Basic(A,newtype,MAT_REUSE_MATRIX,&B);
107:   }

109:   if (reuse == MAT_INPLACE_MATRIX) {
110:     MatHeaderReplace(A,&B);
111:   } else *newmat = B;
112:   return(0);
113: }