qdldl.c in scs-0.3.0

- old
+ new

@@ -1,37 +1,11 @@
 #include "qdldl.h"
-#include "ctrlc.h"
 
 #define QDLDL_UNKNOWN (-1)
 #define QDLDL_USED (1)
 #define QDLDL_UNUSED (0)
 
-// //DEBUG
-// #include <stdio.h>
-// void qdprint_arrayi(const QDLDL_int* data, QDLDL_int n,char* varName){
-
-//   QDLDL_int i;
-//   printf("%s = [",varName);
-//   for(i=0; i< n; i++){
-//     printf("%lli,",data[i]);
-//   }
-//   printf("]\n");
-
-// }
-
-// void qdprint_arrayf(const QDLDL_float* data, QDLDL_int n, char* varName){
-
-//   QDLDL_int i;
-//   printf("%s = [",varName);
-//   for(i=0; i< n; i++){
-//     printf("%.3g,",data[i]);
-//   }
-//   printf("]\n");
-
-// }
-// // END DEBUG
-
 /* Compute the elimination tree for a quasidefinite matrix
    in compressed sparse column form.
 */
 
 QDLDL_int QDLDL_etree(const QDLDL_int  n,
@@ -39,11 +13,11 @@
                       const QDLDL_int* Ai,
                       QDLDL_int* work,
                       QDLDL_int* Lnz,
                       QDLDL_int* etree){
 
-  QDLDL_int sumLnz = 0;
+  QDLDL_int sumLnz;
   QDLDL_int i,j,p;
 
 
   for(i = 0; i < n; i++){
   // zero out Lnz and work.  Set all etree values to unknown
@@ -74,12 +48,23 @@
       }
     }
   }
 
   //compute the total nonzeros in L.  This much
-  //space is required to store Li and Lx
-  for(i = 0; i < n; i++){sumLnz += Lnz[i];}
+  //space is required to store Li and Lx.  Return
+  //error code -2 if the nonzero count will overflow
+  //its unteger type.
+  sumLnz  = 0;
+  for(i = 0; i < n; i++){
+    if(sumLnz > QDLDL_INT_MAX - Lnz[i]){
+      sumLnz = -2;
+      break;
+    }
+    else{
+      sumLnz += Lnz[i];
+    }
+  }
 
   return sumLnz;
 }
 
 
@@ -137,14 +122,10 @@
   Dinv[0] = 1/D[0];
 
   //Start from 1 here. The upper LH corner is trivially 0
   //in L b/c we are only computing the subdiagonal elements
   for(k = 1; k < n; k++){
-    if(scs_is_interrupted()) {
-      scs_printf("interrupt detected in factorization\n");
-      return -1;
-    }
 
     //NB : For each k, we compute a solution to
     //y = L(0:(k-1),0:k-1))\b, where b is the kth
     //column of A that sits above the diagonal.
     //The solution y is then the kth row of L,
@@ -256,30 +237,33 @@
                   const QDLDL_int*   Lp,
                   const QDLDL_int*   Li,
                   const QDLDL_float* Lx,
                   QDLDL_float* x){
 
-QDLDL_int i,j;
+  QDLDL_int i,j;
   for(i = 0; i < n; i++){
-      for(j = Lp[i]; j < Lp[i+1]; j++){
-          x[Li[j]] -= Lx[j]*x[i];
-      }
+    QDLDL_float val = x[i];
+    for(j = Lp[i]; j < Lp[i+1]; j++){
+      x[Li[j]] -= Lx[j]*val;
+    }
   }
 }
 
 // Solves (L+I)'x = b
 void QDLDL_Ltsolve(const QDLDL_int    n,
                    const QDLDL_int*   Lp,
                    const QDLDL_int*   Li,
                    const QDLDL_float* Lx,
                    QDLDL_float* x){
 
-QDLDL_int i,j;
+  QDLDL_int i,j;
   for(i = n-1; i>=0; i--){
-      for(j = Lp[i]; j < Lp[i+1]; j++){
-          x[i] -= Lx[j]*x[Li[j]];
-      }
+    QDLDL_float val = x[i];
+    for(j = Lp[i]; j < Lp[i+1]; j++){
+      val -= Lx[j]*x[Li[j]];
+    }
+    x[i] = val;
   }
 }
 
 // Solves Ax = b where A has given LDL factors
 void QDLDL_solve(const QDLDL_int       n,
@@ -287,12 +271,11 @@
                     const QDLDL_int*   Li,
                     const QDLDL_float* Lx,
                     const QDLDL_float* Dinv,
                     QDLDL_float* x){
 
-QDLDL_int i;
+  QDLDL_int i;
 
-QDLDL_Lsolve(n,Lp,Li,Lx,x);
-for(i = 0; i < n; i++) x[i] *= Dinv[i];
-QDLDL_Ltsolve(n,Lp,Li,Lx,x);
-
+  QDLDL_Lsolve(n,Lp,Li,Lx,x);
+  for(i = 0; i < n; i++) x[i] *= Dinv[i];
+  QDLDL_Ltsolve(n,Lp,Li,Lx,x);
 }