factorns.c

来自「电力系统分析计算 学习调试程序 UNIX / LINUX / CYGWIN 系统」· C语言 代码 · 共 616 行 · 第 1/2 页

        INDEX GetLowval()
        {
          INDEX Iget;
          VALENCE ValGet;
          /* BEGIN */
          ValGet = 0;
          do {
            Iget = ValClass[ValGet];
            ValGet++;
          } while (Iget == 0);
          UNLinkVal(Iget);
          return(Iget);
        }

#ifdef ANSIPROTO
          VALENCE ColumnValence(INDEX Col)
#else
          VALENCE ColumnValence(Col)
          INDEX Col;
#endif
          {
            VALENCE TempCount;
            SparseMatrixElement *ColValPtr;
            /* BEGIN */
            TempCount = 0;
            ColValPtr = Matrix->ColHead[Col];
            while (ColValPtr != NULL) {
              TempCount++;
              ColValPtr = ColValPtr->ColNext;
            }
            return(TempCount);
          }

#ifdef ANSIPROTO
        INDEX GetHighest(INDEX Ipivot)
#else
        INDEX GetHighest(Ipivot)
        INDEX Ipivot;
#endif
        {
          INDEX LargestI,SmallestJ;
          VALENCE ColVal;
          SparseMatrixElement *GetPtr;
          ELEMENTVALUETYPE LargestV;
          float AlphaV;
          /* BEGIN */
          LargestV = 0;  LargestI = 0;
          /* Get the largest pivot on row: */
          GetPtr = Matrix->RowHead[Ipivot];
          while ((GetPtr != NULL) && (GetPtr->Col <= PartitionCol[CurrentColPart])) {
            if (ColNew->p[GetPtr->Col] == 0) {
                if (LargestV <= fabs(GetPtr->Value)) {
                  LargestV = fabs(GetPtr->Value);
                  PivotV = LargestV;
                  DiagVal = GetPtr->Value;
                  LargestI = GetPtr->Col;
                }
            }
            GetPtr = GetPtr->RowNext;
          }
          if (Alpha == 1.0) return(LargestI); /* The most robust case */
          SmallestJ = maxN + 1;
          AlphaV = LargestV * Alpha;
          GetPtr = Matrix->RowHead[Ipivot];
          while ((GetPtr != NULL) && (GetPtr->Col <= PartitionCol[CurrentColPart])) {
            if (ColNew->p[GetPtr->Col] == 0) {
              ColVal = ColumnValence(GetPtr->Col);
              if ((ColVal <= SmallestJ) && (fabs(GetPtr->Value) >= AlphaV)) {
                SmallestJ = ColVal;
                PivotV = fabs(GetPtr->Value);
                DiagVal = GetPtr->Value;
                LargestI = GetPtr->Col;
              }
            }
            GetPtr = GetPtr->RowNext;
          }
          return(LargestI);
        } /* END GetHighest */

#ifdef ANSIPROTO
        void UpdateVals(INDEX Jpivot)
#else
        void UpdateVals(Jpivot)
        INDEX Jpivot;
#endif
        {
          SparseMatrixElement * UpdP1;
          INDEX Jord;
          /* BEGIN */
          UpdP1 = Matrix->ColHead[Jpivot];
          while (UpdP1 != NULL) {
            Jord = UpdP1->Row;
            if ((RowNew->p[Jord] == 0) && (Jord >= Ibeg) && (Jord <= Iend)) {
              UNLinkVal(Jord);
              FindVal(Jord);
              LinkVal(Jord,Valence[Jord]);
            } /* END IF */
            UpdP1 = UpdP1->ColNext;
          } /* END WHILE */
        } /* END UpdateVals */

#ifdef ANSIPROTO
      int OrderMatrix(INDEX Ibeg,INDEX Iend)
#else
      int OrderMatrix(Ibeg,Iend)
      INDEX Ibeg,Iend;
#endif
      {
        int LoopLimit;
        INDEX Istop,Ipivot,Jpivot;
        /* BEGIN OrderMatrix */
        FormClass(Ibeg,Iend);
        Istop = Iend;
        if (Istop > Matrix->n1) Istop = Matrix->n1;
        if (Istop > Matrix->n2) Istop = Matrix->n2;
        if (Istop > Nstop) Istop = Nstop;
        while (Iorder < Istop) {
          LoopLimit = Matrix->n1;
          Iorder++;

          if (Iorder > PartitionCol[CurrentColPart]) CurrentColPart++;

          do {
            LoopLimit--;
            Ipivot = GetLowval();
            Jpivot = GetHighest(Ipivot);
            if (Jpivot == 0) {
              UNLinkVal(Ipivot);
              if (Matrix->n1 == Matrix->n2) {
                fCustomPrint(stderr,"\nUnable to find a nonzero pivot for row %d\n",RowOld->p[Ipivot]);
                fCustomPrint(stderr,"Matrix is probably numerically singular\n");
                return(WARNINGEXIT);
              }
              Valence[Ipivot] = Valence[Ipivot] + DeltaValence;
              LinkVal(Ipivot,Valence[Ipivot]);
            }
          } while (((Jpivot <= 0) || (PivotV <= SINGULARITYZERO)) &&
                  (LinkCount > 0) && (LoopLimit >= 0));

          if (Jpivot <= 0) {
            fCustomPrint(stderr,"\nUnable to find an available pivot for row %d\n",RowOld->p[Ipivot]);
            fCustomPrint(stderr,"Matrix is probably topologically singular\n");
            return(WARNINGEXIT);
            /* Iorder = Istop; */
          } else {
            RowNew->p[Ipivot] = Iorder;
            RowOld->p[Iorder] = Ipivot;
            if (NearZero(PivotV)) {
              fCustomPrint(stderr,"\nEquation %d Depends on Other Equations\n",RowOld->p[Ipivot]);
              fCustomPrint(stderr,"Matrix is probably numerically singular\n");
              return(WARNINGEXIT);
              /* Iorder = Istop; */
            } else {
              ColNew->p[Jpivot] = Iorder;
              ColOld->p[Iorder] = Jpivot;
              SimElim(Ipivot,Jpivot,&DiagVal);
              UpdateVals(Jpivot);
            }
          }
        } /* END WHILE */
        return(0);
      } /* END OrderMatrix */

/* ============================ InvertNormalize ======================== */
  void InvertNormalize()
  {
    INDEX I,Inew;
    SparseMatrixElement *Ptr1;
    ELEMENTVALUETYPE DiagValue;
    /* BEGIN {InvertNormalize} */
    for (Inew=1; Inew<=Nstop; Inew++) {
      I = RowOld->p[Inew];
      Ptr1 = Matrix->RowHead[I];
      while (Ptr1 != NULL) {
        if (RowNew->p[I] == ColNew->p[Ptr1->Col]) {
          /* InvertValue(Ptr1->Value); */
          /* EquateValues(DiagValue,Ptr1->Value); */
          Ptr1->Value = 1.0 / Ptr1->Value;
          DiagValue = Ptr1->Value;
        }
        Ptr1 = Ptr1->RowNext;
      }
      Ptr1 = Matrix->RowHead[I];
      while (Ptr1 != NULL) {
        if (ColNew->p[Ptr1->Col] > RowNew->p[I]) {
          /* MultiplyValues(Ptr1^.Value,Ptr1^.Value,DiagValue); */
          Ptr1->Value = Ptr1->Value * DiagValue;
          Nmult++;
        }
        Ptr1 = Ptr1->RowNext;
      }
    }
  } /* END {InvertNormalize} */

/* ========================== FinishVectors =========================== */
  void FinishVectors()
  {
    INDEX I,K;
    /* BEGIN */
    K = 1;
    while ((RowOld->p[K] != 0) && (K <= Matrix->n1)) K++;
    for (I=1; I<=Matrix->n1; I++) {
      if (RowNew->p[I] == 0) {
        RowNew->p[I] = K;
        K++;
      }
    }
    for (I=1; I<=Matrix->n1; I++) {
      K = RowNew->p[I];
      RowOld->p[K] = I;
    }
    K = 1;
    while ((ColOld->p[K] != 0) && (K <= Matrix->n2)) K++;
    for (I=1; I<=Matrix->n2; I++) {
      if (ColNew->p[I] == 0) {
        ColNew->p[I] = K;
        K++;
      }
    }
    for (I=1; I<=Matrix->n2; I++) ColOld->p[ColNew->p[I]] = I;
  }

/* ============================= factorns ================================ */
#ifdef ANSIPROTO
int factorns(SparseMatrix *Mptr,double Param,IntegerVector *PartRow,IntegerVector *PartCol,
             IntegerVector *P1Row,IntegerVector *P1Col,IntegerVector *P2Row,IntegerVector *P2Col)
#else
int factorns(Mptr,Param,PartRow,PartCol,P1Row,P1Col,P2Row,P2Col)
SparseMatrix *Mptr;
double Param;
IntegerVector *PartRow,*PartCol,*P1Row,*P1Col,*P2Row,*P2Col;
#endif
{ /* FactorNS */
  double MinAlpha = 0.0;
  double MaxAlpha = 1.0;
  SparseMatrixElement *Ptr;
  INDEX i;
  FILE *OutFile;

  /* BEGIN */
  RowNew=P1Row; ColNew=P1Col;
  RowOld=P2Row; ColOld=P2Col;
  Alpha=Param;
  if(Alpha<MinAlpha) Alpha=MinAlpha;
  if(Alpha>MaxAlpha) Alpha=MaxAlpha;
  Matrix = Mptr;
  InitializeLSWR();
  AllocateOrderArrays();
  if (Matrix->n1 < Matrix->n2) Nstop = Matrix->n1;
                          else Nstop = Matrix->n2;
  PartitionAt = PartRow->p;
  Nparts = PartRow->N;
  PartitionCol = PartCol->p;
  NpartCol = PartCol->N;
  CurrentColPart = 1;
  Iorder = 0;
  LinkCount = 0;
  Nfills = 0;
  Nmult = 0;
  for (Ipart=1; Ipart<=Nparts; Ipart++) {
    Ibeg = PartitionAt[Ipart-1]+1;
    Iend = PartitionAt[Ipart];
    if (OrderMatrix(Ibeg,Iend)) return(1);
  }
/*  fCustomPrint(stderr,"  Factorization Fills: %d",Nfills);*/
  FinishVectors();   /* For the case of rectangular matrix factorization */
  InvertNormalize(); /* To convert to true LDU factors, as expected by REPSOL */
/*  fCustomPrint(stderr,"    Multiplications: %d\n",Nmult);*/
  for (i=1; i<=Matrix->n1; i++) {
    Ptr=Matrix->RowHead[i];
    while(Ptr!=NULL){
      Ptr->Row=RowNew->p[Ptr->Row];
      Ptr->Col=ColNew->p[Ptr->Col];
      Ptr=Ptr->RowNext;
    }
  }
#ifdef WINDOWS
  delete[] Lswr;
  delete[] Swr;
  delete[] Valence;
  delete[] FwdLink;
  delete[] BckLink;
  delete[] ValClass;
#else
  free(Lswr);
  free(Swr);
  free(Valence);
  free(FwdLink);
  free(BckLink);
  free(ValClass);
#endif
/*
  if (ExistParameter('d')) {
    OutFile=(FILE *) OpenOutput("prow.dat");
    fCustomPrint(OutFile,"%d\n",RowOld->N);
    for (i=1; i<=RowOld->N; i++) fCustomPrint(OutFile,"%d ",RowOld->p[i]);
    fCustomPrint(OutFile,"\n");
    fclose(OutFile);
    OutFile=(FILE *) OpenOutput("pcol.dat");
    fCustomPrint(OutFile,"%d\n",ColOld->N);
    for (i=1; i<=ColOld->N; i++) fCustomPrint(OutFile,"%d ",ColOld->p[i]);
    fCustomPrint(OutFile,"\n");
    fclose(OutFile);
  }
*/
  return(0);
} /* FactorNS */