Actual source code: baijfact3.c

      1: #define PETSCMAT_DLL

      3: /*
      4:     Factorization code for BAIJ format.
      5: */
     6:  #include src/mat/impls/baij/seq/baij.h
     7:  #include src/inline/ilu.h

      9: /*
     10:     The symbolic factorization code is identical to that for AIJ format,
     11:   except for very small changes since this is now a SeqBAIJ datastructure.
     12:   NOT good code reuse.
     13: */
     16: PetscErrorCode MatLUFactorSymbolic_SeqBAIJ(Mat A,IS isrow,IS iscol,MatFactorInfo *info,Mat *B)
     17: {
     18:   Mat_SeqBAIJ    *a = (Mat_SeqBAIJ*)A->data,*b;
     19:   IS             isicol;
     21:   PetscInt       *r,*ic,i,n = a->mbs,*ai = a->i,*aj = a->j;
     22:   PetscInt       *ainew,*ajnew,jmax,*fill,*ajtmp,nz,bs = A->bs,bs2=a->bs2;
     23:   PetscInt       *idnew,idx,row,m,fm,nnz,nzi,reallocs = 0,nzbd,*im;
     24:   PetscReal      f = 1.0;

     27:   if (A->M != A->N) SETERRQ(PETSC_ERR_ARG_WRONG,"matrix must be square");
     28:   ISInvertPermutation(iscol,PETSC_DECIDE,&isicol);
     29:   ISGetIndices(isrow,&r);
     30:   ISGetIndices(isicol,&ic);

     32:   f = info->fill;
     33:   /* get new row pointers */
     34:   PetscMalloc((n+1)*sizeof(PetscInt),&ainew);
     35:   ainew[0] = 0;
     36:   /* don't know how many column pointers are needed so estimate */
     37:   jmax     = (PetscInt)(f*ai[n] + 1);
     38:   PetscMalloc((jmax)*sizeof(PetscInt),&ajnew);
     39:   /* fill is a linked list of nonzeros in active row */
     40:   PetscMalloc((2*n+1)*sizeof(PetscInt),&fill);
     41:   im       = fill + n + 1;
     42:   /* idnew is location of diagonal in factor */
     43:   PetscMalloc((n+1)*sizeof(PetscInt),&idnew);
     44:   idnew[0] = 0;

     46:   for (i=0; i<n; i++) {
     47:     /* first copy previous fill into linked list */
     48:     nnz     = nz    = ai[r[i]+1] - ai[r[i]];
     49:     if (!nz) SETERRQ(PETSC_ERR_MAT_LU_ZRPVT,"Empty row in matrix");
     50:     ajtmp   = aj + ai[r[i]];
     51:     fill[n] = n;
     52:     while (nz--) {
     53:       fm  = n;
     54:       idx = ic[*ajtmp++];
     55:       do {
     56:         m  = fm;
     57:         fm = fill[m];
     58:       } while (fm < idx);
     59:       fill[m]   = idx;
     60:       fill[idx] = fm;
     61:     }
     62:     row = fill[n];
     63:     while (row < i) {
     64:       ajtmp = ajnew + idnew[row] + 1;
     65:       nzbd  = 1 + idnew[row] - ainew[row];
     66:       nz    = im[row] - nzbd;
     67:       fm    = row;
     68:       while (nz-- > 0) {
     69:         idx = *ajtmp++;
     70:         nzbd++;
     71:         if (idx == i) im[row] = nzbd;
     72:         do {
     73:           m  = fm;
     74:           fm = fill[m];
     75:         } while (fm < idx);
     76:         if (fm != idx) {
     77:           fill[m]   = idx;
     78:           fill[idx] = fm;
     79:           fm        = idx;
     80:           nnz++;
     81:         }
     82:       }
     83:       row = fill[row];
     84:     }
     85:     /* copy new filled row into permanent storage */
     86:     ainew[i+1] = ainew[i] + nnz;
     87:     if (ainew[i+1] > jmax) {

     89:       /* estimate how much additional space we will need */
     90:       /* use the strategy suggested by David Hysom <hysom@perch-t.icase.edu> */
     91:       /* just double the memory each time */
     92:       PetscInt maxadd = jmax;
     93:       /* maxadd = (int)((f*(ai[n]+1)*(n-i+5))/n); */
     94:       if (maxadd < nnz) maxadd = (n-i)*(nnz+1);
     95:       jmax += maxadd;

     97:       /* allocate a longer ajnew */
     98:       PetscMalloc(jmax*sizeof(PetscInt),&ajtmp);
     99:       PetscMemcpy(ajtmp,ajnew,ainew[i]*sizeof(PetscInt));
    100:       PetscFree(ajnew);
    101:       ajnew = ajtmp;
    102:       reallocs++; /* count how many times we realloc */
    103:     }
    104:     ajtmp = ajnew + ainew[i];
    105:     fm    = fill[n];
    106:     nzi   = 0;
    107:     im[i] = nnz;
    108:     while (nnz--) {
    109:       if (fm < i) nzi++;
    110:       *ajtmp++ = fm;
    111:       fm       = fill[fm];
    112:     }
    113:     idnew[i] = ainew[i] + nzi;
    114:   }

    116: #if defined(PETSC_USE_DEBUG)
    117:   if (ai[n] != 0) {
    118:     PetscReal af = ((PetscReal)ainew[n])/((PetscReal)ai[n]);
    119:     PetscLogInfo((A,"MatLUFactorSymbolic_SeqBAIJ:Reallocs %D Fill ratio:given %g needed %g\n",reallocs,f,af));
    120:     PetscLogInfo((A,"MatLUFactorSymbolic_SeqBAIJ:Run with -pc_lu_fill %g or use \n",af));
    121:     PetscLogInfo((A,"MatLUFactorSymbolic_SeqBAIJ:PCLUSetFill(pc,%g);\n",af));
    122:     PetscLogInfo((A,"MatLUFactorSymbolic_SeqBAIJ:for best performance.\n"));
    123:   } else {
    124:     PetscLogInfo((A,"MatLUFactorSymbolic_SeqBAIJ:Empty matrix.\n"));
    125:   }
    126: #endif

    128:   ISRestoreIndices(isrow,&r);
    129:   ISRestoreIndices(isicol,&ic);

    131:   PetscFree(fill);

    133:   /* put together the new matrix */
    134:   MatCreate(A->comm,B);
    135:   MatSetSizes(*B,bs*n,bs*n,bs*n,bs*n);
    136:   MatSetType(*B,A->type_name);
    137:   MatSeqBAIJSetPreallocation_SeqBAIJ(*B,bs,MAT_SKIP_ALLOCATION,PETSC_NULL);
    138:   PetscLogObjectParent(*B,isicol);
    139:   b = (Mat_SeqBAIJ*)(*B)->data;
    140:   b->singlemalloc = PETSC_FALSE;
    141:   PetscMalloc((ainew[n]+1)*sizeof(MatScalar)*bs2,&b->a);
    142:   b->j          = ajnew;
    143:   b->i          = ainew;
    144:   b->diag       = idnew;
    145:   b->ilen       = 0;
    146:   b->imax       = 0;
    147:   b->row        = isrow;
    148:   b->col        = iscol;
    149:   b->pivotinblocks = (info->pivotinblocks) ? PETSC_TRUE : PETSC_FALSE;
    150:   PetscObjectReference((PetscObject)isrow);
    151:   PetscObjectReference((PetscObject)iscol);
    152:   b->icol       = isicol;
    153:   PetscMalloc((bs*n+bs)*sizeof(PetscScalar),&b->solve_work);
    154:   /* In b structure:  Free imax, ilen, old a, old j.
    155:      Allocate idnew, solve_work, new a, new j */
    156:   PetscLogObjectMemory(*B,(ainew[n]-n)*(sizeof(PetscInt)+sizeof(MatScalar)));
    157:   b->maxnz = b->nz = ainew[n];
    158:
    159:   (*B)->factor                 = FACTOR_LU;
    160:   (*B)->info.factor_mallocs    = reallocs;
    161:   (*B)->info.fill_ratio_given  = f;
    162:   if (ai[n] != 0) {
    163:     (*B)->info.fill_ratio_needed = ((PetscReal)ainew[n])/((PetscReal)ai[n]);
    164:   } else {
    165:     (*B)->info.fill_ratio_needed = 0.0;
    166:   }
    167:   return(0);
    168: }