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     ! 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