bug 2003 fixed
Serge Steer [Mon, 7 Jul 2008 16:50:53 +0000 (16:50 +0000)]
scilab/modules/sparse/sci_gateway/fortran/sci_ludel.f
scilab/modules/sparse/sci_gateway/fortran/sci_lufact.f
scilab/modules/sparse/sci_gateway/fortran/sci_luget.f
scilab/modules/sparse/sci_gateway/fortran/sci_lusolve.f
scilab/modules/sparse/src/c/lu.c
scilab/modules/sparse/src/c/lu.h
scilab/modules/sparse/tests/nonreg_tests/bug_2003.tst [new file with mode: 0644]

index e2ae499..b060e94 100644 (file)
@@ -8,7 +8,7 @@ c are also available at
 c http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
       subroutine intludel(id)
       include 'stack.h'
-      integer id(nsiz),top0
+      integer id(nsiz),top0,hand
       integer iadr, sadr
 
       iadr(l)=l+l-1
@@ -37,8 +37,14 @@ c
          return
       endif
       l1 = sadr(il1+4)
+      hand=stk(l1)
+      call ludel1(hand,ierr)
+      if (ierr.ne.0) then
+         err=1
+         call error(247)
+         return
+      endif
 
-      call ludel1(stk(l1))
       if (err .gt. 0) return
 c     
       top=top-rhs
index 1b52426..96504dc 100644 (file)
@@ -10,7 +10,7 @@ c http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
       subroutine intlufact(id)
       include 'stack.h'
       double precision abstol,reltol
-      integer id(nsiz),top0,tops,rhss
+      integer id(nsiz),top0,tops,rhss,fmatindex
       integer iadr, sadr
 c
       iadr(l)=l+l-1
@@ -82,17 +82,10 @@ c
          return
       endif
 
-      lw5=lw
-      lw=lw+1
-      err=lw-lstk(bot)
-      if (err .gt. 0) then
-         call error(17)
-         return
-      endif
-c     
+     
       mx=max(m,n)
       call lufact1(stk(l),istk(il1+5),istk(il1+5+m),mx,nel,
-     $     stk(lw5),abstol,reltol,nrank,ierr)
+     $     fmatindex,abstol,reltol,nrank,ierr)
       if(ierr.gt.0) then
          buf='not enough memory'
          call error(9999)
@@ -110,7 +103,7 @@ c
       istk(il+2)=n
       istk(il+3)=it
       l=sadr(il+4)
-      stk(l)=stk(lw5)
+      stk(l)=fmatindex
       lstk(top+1)=l+1
 c     
       if(lhs .eq.2) then
index d98ea72..4e1ce7b 100644 (file)
@@ -8,8 +8,7 @@ c are also available at
 c http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
       subroutine intluget(id)
       include 'stack.h'
-      double precision ptr
-      integer id(nsiz),top0
+      integer id(nsiz),top0,ptr
       integer iadr, sadr
 c
       iadr(l)=l+l-1
@@ -41,7 +40,12 @@ c
       it1=istk(il1+3)
       l1 = sadr(il1+4)
       ptr=stk(l1)
-      call lusiz1(ptr,nl,nu)
+      call lusiz1(ptr,nl,nu,ierr)
+      if (ierr.ne.0) then
+         err=1
+         call error(247)
+         return
+      endif
       ilp=il1
       lp=sadr(ilp+5+m+m)
       lw=lp+m*(it1+1)
@@ -95,8 +99,12 @@ c
       istk(ilq+4)=n
 c     
       call luget1(ptr,istk(ilp+5),stk(lp),istk(ill+5),stk(ll),
-     $     istk(ilu+5),stk(lu),istk(ilq+5),stk(lq))
-      
+     $     istk(ilu+5),stk(lu),istk(ilq+5),stk(lq),ierr)
+      if (ierr.ne.0) then
+         err=1
+         call error(247)
+         return
+      endif
       return
       end
 c                      ======================================
index cd0d66a..12370e6 100644 (file)
@@ -9,8 +9,8 @@ c http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
 
       subroutine intlusolve(id)
       include 'stack.h'
-      double precision abstol,reltol,hand
-      integer id(nsiz),top0
+      double precision abstol,reltol
+      integer id(nsiz),top0,hand
       integer iadr, sadr
       logical fact
 
@@ -105,14 +105,19 @@ c     b is full
          endif
 c        
          do 40 j=0,n2-1
-            call lusolve1(hand,stk(l2+j*m2),stk(lw3+j*m2))
+            call lusolve1(hand,stk(l2+j*m2),stk(lw3+j*m2),ierr)
             if(it2.eq.1) then
-               call lusolve1(hand,stk(l2i+j*m2),stk(lw3i+j*m2))
+               call lusolve1(hand,stk(l2i+j*m2),stk(lw3i+j*m2),ierr)
+            endif
+            if (ierr.ne.0) then
+               err=1
+               call error(247)
+               return
             endif
             if (err .gt. 0) return
  40      continue
 
-         if(.not.fact) call ludel1(hand)
+         if(.not.fact) call ludel1(hand,ierr)
 c     
          top=top-rhs
          lw0=lw
index d002000..7038494 100644 (file)
@@ -12,6 +12,7 @@
  * the pointer transmitted to f_ is an istk(il1) it can in fact contain
  * something as long as a double
  * Copyright ENPC (Chancelier)
  */
 
 
@@ -60,17 +61,22 @@ static void spGetNumRank(char* eMatrix, int *n)
 }
 
 void C2F(lufact1)(double *val, int *lln, int *col, int *n, int *nel,
-                 long *fmat, double *eps, double *releps, int *nrank, int *ierr)
+                 int *fmatindex, double *eps, double *releps, int *nrank, int *ierr)
 {
   int error,i,i0,i1,k,j;
+  char *fmat;
   spREAL *pelement;
   *ierr = 0;
-  *fmat = (long)spCreate(*n,0,&error);
+  fmat = spCreate(*n,0,&error);
   if (error != spOKAY) {
     *ierr = 1;
     return;
   }
-
+  *fmatindex=addluptr (fmat);
+  if ( *fmatindex == -1) {
+    *ierr = 1;
+    return;
+  }
 
   i0=0;
   i1=i0;
@@ -83,7 +89,9 @@ void C2F(lufact1)(double *val, int *lln, int *col, int *n, int *nel,
       i0=i0+1;
     }
     j=col[k];
-    pelement = spGetElement((char*) *fmat,i,j);
+   
+    pelement = spGetElement(fmat,i,j);
     if (pelement == 0) {
       *ierr=2;
       return;
@@ -92,10 +100,13 @@ void C2F(lufact1)(double *val, int *lln, int *col, int *n, int *nel,
 
   }
   /* Fix the AbsThresold with scilex %eps */
-  spFixThresold((char*) *fmat,*eps,*releps);
+  spFixThresold(fmat,*eps,*releps);
+
   /* spPrint((char *) *fmat,1,1,1); */
-  error = spFactor((char*) *fmat);
-  spGetNumRank((char *) *fmat,nrank);
+  error = spFactor(fmat);
+
+  spGetNumRank(fmat,nrank);
+
   switch (error) {
   case spZERO_DIAG:
     cerro(_("zero_diag: A zero was encountered on the diagonal the matrix "));
@@ -110,6 +121,7 @@ void C2F(lufact1)(double *val, int *lln, int *col, int *n, int *nel,
     *ierr=-2; /* matrix is singular at precision level */
     break;
   }
+
 }
 
 /*
@@ -118,9 +130,15 @@ void C2F(lufact1)(double *val, int *lln, int *col, int *n, int *nel,
  *   b,v
  *      two arrays of size n the matrix size
  */
-void C2F(lusolve1)(long *fmat, double *b, double *x)
+void C2F(lusolve1)(int *fmatindex, double *b, double *x, int *ierr)
 {
-  spSolve((char*) *fmat,(spREAL*)b,(spREAL*)x);
+  char *fmat;
+  if (getluptr((int)*fmatindex, &fmat)==-1){
+    *ierr=1;
+    return;
+  }
+  *ierr = 0;
+  spSolve(fmat,(spREAL*)b,(spREAL*)x);
 }
 
 /*
@@ -128,9 +146,17 @@ void C2F(lusolve1)(long *fmat, double *b, double *x)
  *   *fmat : a pointer to the sparse matrix factored by lufact
  */
 
-void C2F(ludel1)(long *fmat)
+void C2F(ludel1)(int *fmatindex, int *ierr)
 {
-  spDestroy((char*) *fmat);
+  char *fmat;
+  if (getluptr((int)*fmatindex, &fmat)==-1){
+    *ierr=1;
+    return;
+  }
+  *ierr = 0;
+  removeluptr ((int)*fmatindex);
+  spDestroy(fmat);
 }
 
 /*
@@ -286,11 +312,17 @@ static void spLuget(char *eMatrix, int *indP, double *P, int* indl,
 
 
 
-void C2F(luget1)(long *fmat, int *indP, double *P,
+void C2F(luget1)(int *fmatindex, int *indP, double *P,
                 int *indl, double *l, int *indu, double *u,
-                int *indQ, double *Q)
+                int *indQ, double *Q, int *ierr)
 {
-  spLuget((char *) *fmat,indP,P,indl,l,indu,u,indQ,Q);
+  char *fmat;
+  if (getluptr((int)*fmatindex, &fmat)==-1){
+    *ierr=1;
+    return;
+  }
+  *ierr = 0;
+  spLuget(fmat,indP,P,indl,l,indu,u,indQ,Q);
 }
 
 
@@ -323,7 +355,87 @@ static void spLusiz(char *eMatrix, int *lsize,int *usize)
     };
 }
 
-void C2F(lusiz1)(long* fmat, int* lsize, int* usize)
+void C2F(lusiz1)(int* fmatindex, int* lsize, int* usize, int *ierr)
+{
+  char *fmat;
+  if (getluptr((int)*fmatindex, &fmat)==-1){
+    *ierr = 1;
+    return;
+  }  
+  *ierr = 0;
+  spLusiz(fmat,lsize,usize);
+}
+
+char **sci_luptr_table = NULL;
+int sci_luptr_table_size = 0;/* allocated size for pointer table*/
+int sci_luptr_index = 0;/* max index used (one based)*/
+
+
+/**addluptr
+ * This function adds a pointer on a sparse lu factorization to Scilab internal table
+ */
+int addluptr (char *ptr)
+{
+  int i,sel;
+  int rsize=10;
+  if (sci_luptr_table_size==0) { /* first call alloacte a small array of pointers*/
+    sci_luptr_table = (char **)MALLOC(rsize*sizeof(char *));
+    if (sci_luptr_table==NULL) return -1;
+    sci_luptr_table_size += 10;
+  }
+  /* look for a free cell in sci_luptr_table*/
+  sel = -1;
+  for (i=0; i<sci_luptr_index; i++) {
+    if ( sci_luptr_table[i]==NULL) {
+      sel = i;
+      break;
+    }
+  }
+  if (sel == -1) {
+    if (sci_luptr_index < sci_luptr_table_size) {
+      sel=sci_luptr_index++;
+    }
+    else {
+      sci_luptr_table=(char **)REALLOC(sci_luptr_table,(sci_luptr_table_size+rsize)*sizeof(char *));
+      if (sci_luptr_table==NULL) return -1;
+      sci_luptr_table_size += 10;
+      sel=sci_luptr_index++;
+    }
+  }
+  sci_luptr_table[sel] = ptr;
+  return sel+1;
+}
+/**getluptr 
+ * this function returns a pointer on a sparse lu factorization 
+ * given its index (one based) in the table
+ */
+int getluptr(int sel, char **ptr)
+{
+  if (sel > sci_luptr_index || sel < 1) return -1;
+  if (sci_luptr_table[sel-1]==NULL) return -1;
+  *ptr=sci_luptr_table[sel-1];
+  return 0;
+}
+
+/**removeluptr 
+ * This function removes a pointer on a sparse lu factorization 
+ * out of  Scilab internal table given its index in the table
+ */
+int removeluptr (int sel)
+{
+  if (sel > sci_luptr_index || sel < 1) return -1;
+  sci_luptr_table[sel-1]=NULL;
+  if (sel == sci_luptr_index) sci_luptr_index--;
+  return 0;
+}
+
+/**resetluptr 
+ * This function reinitialize the Scilab sparse lu pointer table
+ */
+void resetluptr()
 {
-  spLusiz((char *) *fmat,lsize,usize);
+  FREE(sci_luptr_table);
+  **sci_luptr_table = NULL;
+  sci_luptr_table_size = 0;/* allocated size for pointer table*/
+  sci_luptr_index = 0;/* max index used (one based)*/
 }
index a17a267..b216b93 100644 (file)
 #include "cerro.h"
 #include "machine.h"
 #include "localization.h"
+#include "MALLOC.h"
+
+int addluptr (char *ptr); /* */
+int getluptr (int sel, char **ptr);
+int removeluptr (int sel);
+void resetluptr (void); /* to be used to free the lu pointer table */
 
 void C2F(lufact1)(double *val, int *lln, int *col, int *n, int *nel,
-                 long *fmat, double *eps, double *releps, int *nrank, int *ierr);
+                 int *fmatindex, double *eps, double *releps, int *nrank, int *ierr);
 
-void C2F(lusolve1)(long *fmat, double *b, double *x);
+void C2F(lusolve1)(int *fmatindex, double *b, double *x, int *ierr);
 
-void C2F(ludel1)(long *fmat);
+void C2F(ludel1)(int *fmatindex, int *ierr);
 
-void C2F(lusiz1)(long* fmat, int* lsize, int* usize);
+void C2F(lusiz1)(int *fmatindex, int* lsize, int* usize, int *ierr);
 
-void C2F(luget1)(long *fmat, int *indP, double *P,
+void C2F(luget1)(int *fmatindex, int *indP, double *P,
                 int *indl, double *l, int *indu, double *u,
-                int *indQ, double *Q);
+                int *indQ, double *Q, int *ierr);
 #endif /* !__LU_H__ */
diff --git a/scilab/modules/sparse/tests/nonreg_tests/bug_2003.tst b/scilab/modules/sparse/tests/nonreg_tests/bug_2003.tst
new file mode 100644 (file)
index 0000000..d9f5c37
--- /dev/null
@@ -0,0 +1,39 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2007-2008 - INRIA - Serge STEER <serge.steer@inria.fr>
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+
+// <-- Non-regression test for bug 2003 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=2003
+//
+// <-- Short Description -->
+//lusolve crashes scilab when handle has been freed
+
+a=[0.2,0.6,0.6,0.2,0.3;
+0.8,0.8,0.7,0.2,0.9;
+0,0.7,0.7,0.2,0.2;
+0.3,0.9,0.2,0.9,0.3;
+0.7,0.1,0.5,0.7,0.4];
+b=[0.3;0.6;0.5;0.3;0.6];
+A=sparse(a);
+[h,rk]=lufact(A);
+x=lusolve(h,b);
+if norm(a*x-b)>1d-10 then pause,end
+[P,L,U,Q]=luget(h);
+if norm(P*L*U*Q-A)>1d-10 then pause,end
+ludel(h)
+if execstr('x=lusolve(h,b);','errcatch')<>247 then pause,end
+if execstr('ludel(h);','errcatch')<>247 then pause,end
+if execstr('[P,L,U,Q]=luget(h);','errcatch')<>247 then pause,end
+
+//try to allocate a lot of handles
+for k=1:20
+  [h1,rk]=lufact(A);
+end    
+
+x=lusolve(h1,b);
+if norm(a*x-b)>1d-10 then pause,end