[BACK]Return to numerical.c CVS log [TXT][DIR] Up to [local] / OpenXM_contrib2 / asir2000 / builtin

Diff for /OpenXM_contrib2/asir2000/builtin/numerical.c between version 1.3 and 1.7

version 1.3, 2000/08/22 05:03:59 version 1.7, 2018/03/29 01:32:50
Line 45 
Line 45 
  * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,   * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
  * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.   * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
  *   *
  * $OpenXM: OpenXM_contrib2/asir2000/builtin/numerical.c,v 1.2 2000/08/21 08:31:20 noro Exp $   * $OpenXM: OpenXM_contrib2/asir2000/builtin/numerical.c,v 1.6 2003/03/07 06:39:55 noro Exp $
 */  */
 #include "ca.h"  #include "ca.h"
 #include "parse.h"  #include "parse.h"
   
 #if LAPACK  #if defined(LAPACK)
 #include <f2c.h>  #include <f2c.h>
   
 void Pflinsolve();  void Pflinsolve();
   
 struct ftab numerical_tab[] = {  struct ftab numerical_tab[] = {
         {"flinsolve",Pflinsolve,2},    {"flinsolve",Pflinsolve,2},
         {0,0,0},    {0,0,0},
 };  };
   
 void Pflinsolve(arg,rp)  void Pflinsolve(arg,rp)
 NODE arg;  NODE arg;
 VECT *rp;  VECT *rp;
 {  {
         MAT mat;    MAT mat;
         VECT b,r;    VECT b,r;
         int i,j;    int i,j;
         integer n,nrhs,lda,ldb,info;    integer n,nrhs,lda,ldb,info;
         integer *ipiv;    integer *ipiv;
         doublereal *m_array,*b_array;    doublereal *m_array,*b_array;
         Num **mb;    Num **mb;
         Num *vb;    Num *vb;
         Real *rb;    Real *rb;
         char trans;    char trans;
   
         mat=(MAT)ARG0(arg); b=(VECT)ARG1(arg);    mat=(MAT)ARG0(arg); b=(VECT)ARG1(arg);
         n=mat->row;    n=mat->row;
         if ( n != mat->col )    if ( n != mat->col )
                 error("flinsolve : non-square matrix");      error("flinsolve : non-square matrix");
         if ( n != b->len )    if ( n != b->len )
                 error("flinsolve : inconsistent rhs");      error("flinsolve : inconsistent rhs");
         m_array = (doublereal *)MALLOC_ATOMIC(n*n*sizeof(double));    m_array = (doublereal *)MALLOC_ATOMIC(n*n*sizeof(double));
         b_array = (doublereal *)MALLOC_ATOMIC(n*sizeof(double));    b_array = (doublereal *)MALLOC_ATOMIC(n*sizeof(double));
         ipiv = (integer *)MALLOC_ATOMIC(n*sizeof(integer));    ipiv = (integer *)MALLOC_ATOMIC(n*sizeof(integer));
         for ( i = 0, mb = (Num **)mat->body; i < n; i++ )    for ( i = 0, mb = (Num **)mat->body; i < n; i++ )
                 for ( j = 0; j < n; j++ )      for ( j = 0; j < n; j++ )
                         m_array[i*n+j] = ToReal(mb[j][i]);        m_array[i*n+j] = ToReal(mb[j][i]);
         for ( i = 0, vb = (Num *)b->body; i < n; i++ )    for ( i = 0, vb = (Num *)b->body; i < n; i++ )
                         b_array[i] = ToReal(vb[i]);        b_array[i] = ToReal(vb[i]);
         trans = 'N';    trans = 'N';
         nrhs = 1;    nrhs = 1;
         lda = n;    lda = n;
         ldb = n;    ldb = n;
         dgetrf_(&n,&n,m_array,&lda,ipiv,&info);    dgetrf_(&n,&n,m_array,&lda,ipiv,&info);
         if ( info )    if ( info )
                 error("flinsolve : LU factorization failed");      error("flinsolve : LU factorization failed");
         dgetrs_(&trans,&n,&nrhs,m_array,&lda,ipiv,b_array,&ldb,&info);    dgetrs_(&trans,&n,&nrhs,m_array,&lda,ipiv,b_array,&ldb,&info);
         if ( info )    if ( info )
                 error("flinsolve : solving failed");      error("flinsolve : solving failed");
         MKVECT(r,n);    MKVECT(r,n);
         for ( i = 0, rb = (Real *)r->body; i < n; i++ )    for ( i = 0, rb = (Real *)r->body; i < n; i++ )
                 MKReal(b_array[i],rb[i]);      MKReal(b_array[i],rb[i]);
         *rp = r;    *rp = r;
 }  }
 #endif  #endif
Legend:
Removed from v.1.3  
changed lines
  Added in v.1.7
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>