Return to numerical.c CVS log

Up to [local] / OpenXM_contrib2 / asir2000 / builtin

Annotation of OpenXM_contrib2/asir2000/builtin/numerical.c, Revision 1.1.1.1

1.1       noro        1: /* $OpenXM: OpenXM/src/asir99/builtin/numerical.c,v 1.1.1.1 1999/11/10 08:12:25 noro Exp $ */
                      2: #include "ca.h"
                      3: #include "parse.h"
                      4:
                      5: #if LAPACK
                      6: #include <f2c.h>
                      7:
                      8: void Pflinsolve();
                      9:
                     10: struct ftab numerical_tab[] = {
                     11:        {"flinsolve",Pflinsolve,2},
                     12:        {0,0,0},
                     13: };
                     14:
                     15: void Pflinsolve(arg,rp)
                     16: NODE arg;
                     17: VECT *rp;
                     18: {
                     19:        MAT mat;
                     20:        VECT b,r;
                     21:        int i,j;
                     22:        integer n,nrhs,lda,ldb,info;
                     23:        integer *ipiv;
                     24:        doublereal *m_array,*b_array;
                     25:        Num **mb;
                     26:        Num *vb;
                     27:        Real *rb;
                     28:        char trans;
                     29:
                     30:        mat=(MAT)ARG0(arg); b=(VECT)ARG1(arg);
                     31:        n=mat->row;
                     32:        if ( n != mat->col )
                     33:                error("flinsolve : non-square matrix");
                     34:        if ( n != b->len )
                     35:                error("flinsolve : inconsistent rhs");
                     36:        m_array = (doublereal *)MALLOC_ATOMIC(n*n*sizeof(double));
                     37:        b_array = (doublereal *)MALLOC_ATOMIC(n*sizeof(double));
                     38:        ipiv = (integer *)MALLOC_ATOMIC(n*sizeof(integer));
                     39:        for ( i = 0, mb = (Num **)mat->body; i < n; i++ )
                     40:                for ( j = 0; j < n; j++ )
                     41:                        m_array[i*n+j] = ToReal(mb[j][i]);
                     42:        for ( i = 0, vb = (Num *)b->body; i < n; i++ )
                     43:                        b_array[i] = ToReal(vb[i]);
                     44:        trans = 'N';
                     45:        nrhs = 1;
                     46:        lda = n;
                     47:        ldb = n;
                     48:        dgetrf_(&n,&n,m_array,&lda,ipiv,&info);
                     49:        if ( info )
                     50:                error("flinsolve : LU factorization failed");
                     51:        dgetrs_(&trans,&n,&nrhs,m_array,&lda,ipiv,b_array,&ldb,&info);
                     52:        if ( info )
                     53:                error("flinsolve : solving failed");
                     54:        MKVECT(r,n);
                     55:        for ( i = 0, rb = (Real *)r->body; i < n; i++ )
                     56:                MKReal(b_array[i],rb[i]);
                     57:        *rp = r;
                     58: }
                     59: #endif