Changeset 782fcd8 in git

Timestamp:

Jul 6, 2000, 3:32:25 PM (24 years ago)

Author:

Wilfred Pohl <pohl@…>

Branches:

(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')

Children:

80c3d12cb5c52a12ef67dfc3e65568716ca01a42

Parents:

638c809b49d290a8bcb1aa0d3f1ca8d9890f77f4

Message:

linear solving in bareiss


git-svn-id: file:///usr/local/Singular/svn/trunk@4481 2c84dea3-7e68-4137-9b89-c4e89433aadc

Location:

Singular

Files:

• LIB/ntsolve.lib (modified) (14 diffs)
• ideals.cc (modified) (2 diffs)
• iparith.cc (modified) (2 diffs)
• polys1.cc (modified) (2 diffs)
• sparsmat.cc (modified) (18 diffs)

Singular/LIB/ntsolve.lib

@@ -1 +1 @@
 ///////////////////////////////////////////////////////////////////////////////
-version="$Id: ntsolve.lib,v 1.5 2000-01-12 11:01:39 pohl Exp $";
+version="$Id: ntsolve.lib,v 1.6 2000-07-06 13:32:25 pohl Exp $";
 info="
 LIBRARY: ntsolve.lib ONE REAL SOLUTION OF POLYNOMIAL SYSTEMS (NEWTON ITERATION)
@@ -10 +10 @@
 ///////////////////////////////////////////////////////////////////////////////
 
-proc nt_solve( ideal gls, vector ini, intvec ipar )
+proc nt_solve( ideal gls, ideal ini, intvec ipar )
 "USAGE:   nt_solve(gls,ini,ipar);
          gls: the equations
-         ini: the vector of initial values
+         ini: the ideal of initial values
          ipar: control
            ipar[1] - max. number of iterations
@@ -24 +24 @@
 ASSUME:  gls is a zerodimensional ideal with
          nvars(basering) = size(gls) (> 1)
-RETURN:  vector of one solution (if found)
+RETURN:  ideal of one solution (if found)
          0 (else)
 EXAMPLE: example nt_solve; shows an example
@@ -35 +35 @@
     if (size(ini) != di){
       ERROR("wrong number of initial values");}
+    int prec = system("getPrecDigits"); // precision
 
     int i1,i2,i3;
@@ -40 +41 @@
     int itmax, acc, prot;
     if (i1 < 1){itmax = 100;}else{itmax = ipar[1];}
-    if (i1 < 2){acc = 10;}else{acc = ipar[2];}
+    if (i1 < 2){acc = prec/2;}else{acc = ipar[2];}
     if (i1 < 3){prot = 0;}else{prot = ipar[3];}
-
-    int prec = acc+5; // precision in the working ring
+    if ((acc <= 0)||(acc > prec-1)){acc = prec-1;}
+
     int dpl = di+1;
     string out; // for prot != 0 and more
-    intvec permut; // the permutations in bareiss
     out = "ring rnewton=(real,prec),("+varstr(basering)+"),(c,dp);";
     execute(out);
@@ -57 +57 @@
       setring rn; kill rnewton;
       ERROR("one var not in equation");}
-    vector direction,ini1;
+    list direction;
+    ideal ini1;
     ini1 = imap(rn,ini);
     number dum,y1,y2,y3,genau;
@@ -69 +70 @@
       nt = subst(nt,var(i1),ini1[i1]);}
     // now we have in sub the general structure
-    // and in nt the strukture with subst. vars
+    // and in nt the structure with subst. vars
 
     // compute initial error
@@ -81 +82 @@
 
     // find newton direction
-    list bar = bareiss(nt,1,-1);
-    nt = bar[1];
-    permut = bar[2];
-    kill bar;
-    direction=nt[di][dpl]/leadcoef(nt[di][di])*gen(di);
-    for(i1=di-1;i1>0;i1--){
-      y3 = leadcoef(nt[i1][dpl]);
-      for(i2=di;i2>i1;i2--){
-        y3 = y3-leadcoef(direction[i2])*leadcoef(nt[i1][i2]);}
-      direction = direction+(y3/leadcoef(nt[i1][i1])*gen(i1));}
+    direction=bareiss(nt,1,-1);
 
     // find dumping
-    dum = linesearch(gls1,ini1,direction,permut,y1,dum,genau);
+    dum = linesearch(gls1,ini1,direction[1],y1,dum,genau);
+    if (i3%5 == 0)
+    {
+      if (dum <= 0.000001)
+      {
+        dum = 1.0;
+      }
+    }
     if(prot){out="  dumping = "+string(dum);out;}
 
     // new value
     for(i1=di;i1>0;i1--){
-      ini1=ini1-dum*direction[i1]*gen(permut[i1]);}
+      ini1[i1]=ini1[i1]-dum*direction[1][i1];}
     nt = sub;
     for (i1=di;i1>0;i1--){
@@ -107 +106 @@
   }
 
+    if (y1>y2){
+      "WARNING: no convergence";}
     setring rn;
-    if (y1>y2){
-      kill rnewton;
-      ERROR("no convergence");}
     ini = imap(rnewton,ini1);
     kill rnewton;
@@ -118 +116 @@
 {
     "EXAMPLE:";echo=2;
-    ring rsq = (real,16),(x,y,z,w),(c,lp);
+    ring rsq = (real,40),(x,y,z,w),lp;
     ideal gls =  x2+y2+z2-10, y2+z3+w-8, xy+yz+xz+w5 - 1,w3+y;
-    vector ini = [3.1,2.9,1.1,0.5];
-    intvec ipar = 200,8,1;
-    vector sol = nt_solve(gls,ini,ipar);
+    ideal ini = 3.1,2.9,1.1,0.5;
+    intvec ipar = 200,0,1;
+    ideal sol = nt_solve(gls,ini,ipar);
     sol;
 }
@@ -129 +127 @@
 static proc sqrt (number wr, number wa, number wg)
 {
-  number we=wr*wg;
+  number es,we;
   number wb=wa;
   number wf=wb*wb-wr;
-  while (wf>we)
+  if(wf>0){
+    es=wf;}
+  else{
+    es=-wf;}
+  we=wg*es;
+  while (es>we)
   {
     wf=wf/(wb+wb);
     wb=wb-wf;
     wf=wb*wb-wr;
+    if(wf>0){
+      es=wf;}
+    else{
+      es=-wf;}
   }
   return(wb);
-}
-
-static proc vl2norm (vector H, number wg)
-{
-  number wa,wb;
-  int wi;
-  wa = leadcoef(H[1]);
-  wa = wa*wa;
-  for(wi=size(H);wi>1;wi--)
-  {
-    wb=leadcoef(H[wi]);
-    wa=wa+wb*wb;
-  }
-  return(sqrt(wa,wa,wg));
 }
 
@@ -185 +178 @@
 
 static
-proc linesearch(ideal nl, vector aa, vector bb, intvec pe,
+proc linesearch(ideal nl, ideal aa, ideal bb,
 number z1, number tt, number gg)
 {
   int ii,d;
-  vector cc;
-  ideal jn;
-  number z2,z3,e1;
+  ideal cc,jn;
+  number ss,z2,z3,mm;
+
+  mm=0.000001;
+  ss=tt;
   d=size(nl);
   cc=aa;
-  for(ii=d;ii>0;ii--){cc=cc-tt*bb[ii]*gen(pe[ii]);}
+  for(ii=d;ii>0;ii--){cc[ii]=cc[ii]-ss*bb[ii];}
   jn=nl;
   for(ii=d;ii>0;ii--){jn=subst(jn,var(ii),cc[ii]);}
   z2=il2norm(jn,gg);
   z3=-1;
-  e1=1.0e-6;
   while(z2>=z1)
   {
-    if(tt<e1){return (e1);}
-    tt=0.5*tt;
+    ss=0.5*ss;
+    if(ss<mm){return (mm);}
     cc=aa;
     for(ii=d;ii>0;ii--)
     {
-      cc=cc-tt*bb[ii]*gen(pe[ii]);
+      cc[ii]=cc[ii]-ss*bb[ii];
     }
     jn=nl;
@@ -218 +212 @@
     while(z3<z2)
     {
-      tt=tt+tt;
+      ss=ss+ss;
       cc=aa;
       for(ii=d;ii>0;ii--)
       {
-        cc=cc-tt*bb[ii]*gen(pe[ii]);
+        cc[ii]=cc[ii]-ss*bb[ii];
       }
       jn=nl;
@@ -232 +226 @@
   z2=z2-z1;
   z3=z3-z1;
-  tt=0.25*tt*(z3-4*z2)/(z3-2*z2);
-  if(tt<e1){return (e1);}
-  return(tt);
+  ss=0.25*ss*(z3-4*z2)/(z3-2*z2);
+  if(ss>1.0){return (1.0);}
+  if(ss<mm){return (mm);}
+  return(ss);
 }
 ///////////////////////////////////////////////////////////////////////////////

Singular/ideals.cc

@@ -2 +2 @@
 *  Computer Algebra System SINGULAR     *
 ****************************************/
-/* $Id: ideals.cc,v 1.95 2000-05-24 14:32:37 siebert Exp $ */
+/* $Id: ideals.cc,v 1.96 2000-07-06 13:30:01 pohl Exp $ */
 /*
 * ABSTRACT - all basic methods to manipulate ideals
@@ -314 +314 @@
 }
 
+/*2
+*test if the ideal has only constant polynomials
+*/
+BOOLEAN idIsConstant(ideal id)
+{
+  int k;
+  for (k = IDELEMS(id)-1; k>=0; k--)
+  {
+    if (pIsConstantPoly(id->m[k]) == FALSE)
+      return FALSE;
+  }
+  return TRUE;
+}
 
 /*2

Singular/iparith.cc

@@ -2 +2 @@
 *  Computer Algebra System SINGULAR     *
 ****************************************/
-/* $Id: iparith.cc,v 1.216 2000-06-06 15:46:34 Singular Exp $ */
+/* $Id: iparith.cc,v 1.217 2000-07-06 13:30:02 pohl Exp $ */
 
 /*
@@ -3984 +3984 @@
 static BOOLEAN jjBAREISS3(leftv res, leftv u, leftv v, leftv w)
 {
-  lists l=smCallNewBareiss((ideal)u->Data(),(int)v->Data(),(int)w->Data());
+  lists l;
+  int k=(int)w->Data();
+  if (k>=0)
+  {
+    l=smCallNewBareiss((ideal)u->Data(),(int)v->Data(),(int)w->Data());
+  }
+  else
+  {
+    l=smCallSolv((ideal)u->Data());
+  }
   res->data = (char *)l;
   return FALSE;

Singular/polys1.cc

@@ -2 +2 @@
 *  Computer Algebra System SINGULAR     *
 ****************************************/
-/* $Id: polys1.cc,v 1.37 2000-05-05 09:46:57 Singular Exp $ */
+/* $Id: polys1.cc,v 1.38 2000-07-06 13:30:04 pohl Exp $ */
 
 /*
@@ -42 +42 @@
     }
     if (pGetComp(p)) return FALSE;
+  }
+  return TRUE;
+}
+
+/*2
+*test if the polynom is a constant
+*/
+BOOLEAN   pIsConstantPoly(poly p)
+{
+  while(p!=NULL)
+  {
+    for (int i=pVariables;i;i--)
+    {
+      if (pGetExp(p,i)!=0) return FALSE;
+    }
+    pIter(p);
   }
   return TRUE;

Singular/sparsmat.cc

@@ -2 +2 @@
 *  Computer Algebra System SINGULAR     *
 ****************************************/
-/* $Id: sparsmat.cc,v 1.29 2000-06-21 07:32:20 pohl Exp $ */
+/* $Id: sparsmat.cc,v 1.30 2000-07-06 13:30:04 pohl Exp $ */
 
 /*
@@ -97 +97 @@
   void smFinalMult();
   void smSparseHomog();
-  void smRealPivot();
   void smWeights();
   void smPivot();
@@ -472 +471 @@
     if (act < 2)
     {
+      if (TEST_OPT_PROT) PrintS(".\n");
       this->smFinalMult();
       this->smPivDel();
@@ -532 +532 @@
     if (act < y)
     {
+      if (TEST_OPT_PROT) PrintS(".\n");
       this->smCopToRes();
       return;
@@ -543 +544 @@
 *   - the method will finish for number of unreduced columns < y
 */
-void sparse_mat::smNewBareiss(int x, int y0)
-{
-  int y=y0;
+void sparse_mat::smNewBareiss(int x, int y)
+{
   if ((x > 0) && (x < nrows))
   {
@@ -559 +559 @@
   normalize = this->smCheckNormalize();
   if (normalize) this->smNormalize();
-  if(y0>=0)
-    this->smPivot();
-  else
-    this->smRealPivot();
+  this->smPivot();
   this->smSelectPR();
   this->sm1Elim();
@@ -581 +578 @@
   {
     if (normalize) this->smNormalize();
-    if(y0>=0)
-      this->smNewPivot();
-    else
-      this->smRealPivot();
+    this->smNewPivot();
     this->smSelectPR();
     this->smMultCol();
@@ -597 +591 @@
     if (act <= y)
     {
+      if (TEST_OPT_PROT) PrintS(".\n");
       this->smFinalMult();
       this->smCopToRes();
@@ -605 +600 @@
 
 /* ----------------- pivot method ------------------ */
-
-void sparse_mat::smRealPivot()
-{
-  smpoly a;
-  number nopt,n1;
-  int i, copt, ropt;
-
-  nopt=nInit(0);
-  for (i=act; i; i--)
-  {
-    a = m_act[i];
-    loop
-    {
-      if (a->pos > tored)
-        break;
-      n1=pGetCoeff(a->m);
-      if(nGreaterZero(n1))
-      {
-        if(nGreater(n1,nopt))
-        {
-          nDelete(&nopt);
-          nopt=nCopy(n1);
-          copt=i;
-          ropt=a->pos;
-        }
-      }
-      else
-      {
-        n1=nNeg(n1);
-        if(nGreater(n1,nopt))
-        {
-          nDelete(&nopt);
-          nopt=nCopy(n1);
-          copt=i;
-          ropt=a->pos;
-        }
-        n1=nNeg(n1);
-      }
-      a = a->n;
-      if (a == NULL)
-        break;
-    }
-  }
-  rpiv = ropt;
-  cpiv = copt;
-  nDelete(&nopt);
-  if (cpiv != act)
-  {
-    a = m_act[act];
-    m_act[act] = m_act[cpiv];
-    m_act[cpiv] = a;
-  }
-}
 
 /*
@@ -2413 +2355 @@
   void smRowToCol();
   void smSelectPR();
-  void smRealWeights();
   void smRealPivot();
   void smZeroToredElim();
@@ -2434 +2375 @@
 * uses  Gauss-elimination
 */
-ideal smCallSolv(ideal I)
+lists smCallSolv(ideal I)
 {
   sparse_number_mat *linsolv;
@@ -2440 +2381 @@
   ring origR;
   sip_sring tmpR;
-  ideal rr, res;
-
+  ideal rr, ss;
+  lists res;
+
+  if (idIsConstant(I)==FALSE)
+  {
+    WerrorS("symbol in equation");
+    return NULL;
+  }
   I->rank = idRankFreeModule(I);
   if (smCheckSolv(I)) return NULL;
+  res=(lists)AllocSizeOf(slists);
   rr=smRingCopy(I,&origR,tmpR);
-  res = NULL;
+  ss = NULL;
   linsolv = new sparse_number_mat(rr);
   linsolv->smTriangular();
@@ -2451 +2399 @@
   {
     linsolv->smSolv();
-    res = linsolv->smRes2Ideal();
+    ss = linsolv->smRes2Ideal();
   }
   else
     WerrorS("singular problem for linsolv");
   delete linsolv;
-  if ((origR!=NULL) && (res!=NULL))
+  if ((origR!=NULL) && (ss!=NULL))
   {
     rChangeCurrRing(origR,TRUE);
-    rr = idInit(IDELEMS(res), 1);
-    for (k=0;k<IDELEMS(res);k++)
-      rr->m[k] = prCopyR(res->m[k], &tmpR);
+    rr = idInit(IDELEMS(ss), 1);
+    for (k=0;k<IDELEMS(ss);k++)
+      rr->m[k] = prCopyR(ss->m[k], &tmpR);
     rChangeCurrRing(&tmpR,FALSE);
-    idDelete(&res);
-    res = rr;
+    idDelete(&ss);
+    ss = rr;
   }
   if(origR!=NULL)
     smRingClean(origR,tmpR);
+  res->Init(1);
+  res->m[0].rtyp=IDEAL_CMD;
+  res->m[0].data=(void *)ss;
   return res;
 }
@@ -2539 +2490 @@
     if (sing != 0) return;
   }
+  if (TEST_OPT_PROT) PrintS(".\n");
   piv = m_act[1];
   rpiv = piv->pos;
@@ -2646 +2598 @@
 
 /*
-* prepare smPivot, compute weights for rows and columns
-* and the weight for all points
-*/
-void sparse_number_mat::smRealWeights()
-{
-  int wc;
+* compute pivot
+*/
+void sparse_number_mat::smRealPivot()
+{
   smnumber a;
-  int i;
-
-  for (i=tored; i; i--) wrw[i] = 0; // ???
-  for (i=act; i; i--)
-  {
-    wc = 0;
-    a = m_act[i];
-    loop
-    {
-      wc++;
-      wrw[a->pos]++;
-      a = a->n;
-      if ((a == NULL) || (a->pos > tored))
-        break;
-    }
-    wcl[i] = wc;
-  }
-}
-
-/*
-* compute pivot
-*/
-void sparse_number_mat::smRealPivot()
-{
-  int wopt = 1<<30;
-  int w;
-  smnumber a;
-  number x, xo, xu;
+  number x, xo;
   int i, copt, ropt;
 
-  this->smRealWeights();
-  nNew(&xo);
-  nNew(&xu);
+  xo=nInit(0);
   for (i=act; i; i--)
   {
@@ -2691 +2612 @@
     while ((a!=NULL) && (a->pos<=tored))
     {
-      w = (wcl[i]-1)*(wrw[a->pos]-1);
-      if (w==0) // row or column with only one point
-      {
-        rpiv = a->pos;
-        if (i != act)
-        {
-          a = m_act[act];
-          m_act[act] = m_act[i];
-          m_act[i] = a;
-        }
-        return;
-      }
-      if (w == wopt)
-      {
-        x = a->m;
-        if (nGreater(x,xo) || nGreater(xu,x))
-        {
-          nDelete(&xu);
+      x = a->m;
+      if (nGreaterZero(x))
+      {
+        if (nGreater(x,xo))
+        {
           nDelete(&xo);
           xo = nCopy(x);
-          xu = nCopy(x);
-          if (nGreaterZero(xu)) xu = nNeg(xu);
-          else xo = nNeg(xo);
           copt = i;
           ropt = a->pos;
         }
       }
-      else if (w < wopt)
-      {
-        x = a->m;
-        wopt = w;
-        nDelete(&xu);
-        nDelete(&xo);
-        xo = nCopy(x);
-        xu = nCopy(x);
-        if (nGreaterZero(xu)) xu = nNeg(xu);
-        else xo = nNeg(xo);
-        copt = i;
-        ropt = a->pos;
+      else
+      {
+        xo = nNeg(xo);
+        if (nGreater(xo,x))
+        {
+          nDelete(&xo);
+          xo = nCopy(x);
+          copt = i;
+          ropt = a->pos;
+        }
+        xo = nNeg(xo);
       }
       a = a->n;
@@ -2742 +2646 @@
   }
   nDelete(&xo);
-  nDelete(&xu);
 }
 
@@ -2834 +2737 @@
   int i;
 
+  if (TEST_OPT_PROT)
+  {
+    if ((crd+1)%10)
+      PrintS(".");
+    else
+      PrintS(".\n");
+  }
   a = m_act[act];
   if (a->pos < rpiv)