Changeset d98d1b in git

Timestamp:

Oct 25, 2010, 4:09:29 PM (14 years ago)

Author:

Frank Seelisch <seelisch@…>

Branches:

(u'fieker-DuVal', '117eb8c30fc9e991c4decca4832b1d19036c4c65')(u'spielwiese', 'b4f17ed1d25f93d46dbe29e4b499baecc2fd51bb')

Children:

8730345c41037290c50978c95d29abe8d5606708

Parents:

386a8636dfd438b85470e212f44acfc042909ef7

Message:

bug fix in qr double shift

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

File:

• kernel/linearAlgebra.cc (modified) (1 diff)

kernel/linearAlgebra.cc

@@ -1049 +1049 @@
   MATELEM(c,3,1) = ppMult_qq(MATELEM(H,2,1), MATELEM(H,3,2));
   nDelete(&trace); nDelete(&det);
-  
-  matrix uVec; matrix hMat;
-  tmp1 = hessenbergStep(c, uVec, hMat, tolerance); nDelete(&tmp1);
-  /* now replace H by hMat * H * hMat: */
-  matrix wMat = mpMult(hMat, H); idDelete((ideal*)&H);
-  matrix H1 = mpMult(wMat, hMat);
-  idDelete((ideal*)&wMat); idDelete((ideal*)&hMat);
-  /* now need to re-establish Hessenberg form of H1 and put it in H */
-  hessenberg(H1, wMat, H, tolerance);
-  idDelete((ideal*)&wMat); idDelete((ideal*)&H1);
+
+  /* for applying hessenbergStep, we need to make sure that c[1, 1] is
+     not zero */
+  if ((MATELEM(c,1,1) != NULL) &&
+      ((MATELEM(c,2,1) != NULL) || (MATELEM(c,3,1) != NULL)))  
+  {
+    matrix uVec; matrix hMat;
+    tmp1 = hessenbergStep(c, uVec, hMat, tolerance); nDelete(&tmp1);
+    /* now replace H by hMat * H * hMat: */
+    matrix wMat = mpMult(hMat, H); idDelete((ideal*)&H);
+    matrix H1 = mpMult(wMat, hMat);
+    idDelete((ideal*)&wMat); idDelete((ideal*)&hMat);
+    /* now need to re-establish Hessenberg form of H1 and put it in H */
+    hessenberg(H1, wMat, H, tolerance);
+    idDelete((ideal*)&wMat); idDelete((ideal*)&H1);
+  }
+  else if ((MATELEM(c,1,1) == NULL) && (MATELEM(c,2,1) != NULL))
+  {
+    swapRows(1, 2, H);
+    swapColumns(1, 2, H);
+  }
+  else if ((MATELEM(c,1,1) == NULL) && (MATELEM(c,3,1) != NULL))
+  {
+    swapRows(1, 3, H);
+    swapColumns(1, 3, H);
+  }
+  else
+  { /* c is the zero vector or a multiple of e_1;
+       no hessenbergStep needed */ }
 }