数値計算レポート第5回

G99P057-3
斎藤卓也

提出期限: 2001年8月15日
提出日: 2001年8月15日

1  問題

課題は「数値計算用インタプリタを作れ 」である。

2  問題の解析&考察

今回の課題は、flex, bisonを使用して、数値計算を行うことが出来るインタプリタを作れというものであった。数値計算のホームページを見ると、「前田氏」という人が作った、もうかなり完成されている形の計算機があったのだが、どうもプログラムが複雑である。

他人の作ったプログラムを読むのは、なかなか大変なものであり、夏休みに入り、実家に帰省しており、しかもココではプリンタが故障していて印刷さえ出来ないという状況であったので、とても前田氏の電卓プログラムを全て理解する気にはなれなかった。

そこで、前田氏の電卓を改良すれば、簡単に高機能な電卓が作れそうだが、良く分からないまま改良をして機能追加しても、余り勉強にならないので、大したものではなくても良いから、自分でflex, bisonを使って1から作ってみようと思った。

まず、問題はその開発環境である。今は実家なので、Linuxマシンが無い。ノートパソコンにLinuxを入れるという手もあるにはあるが、20GBのHDDをバックアップを取る手段さえ無いため、事実上不可能である。そこで仕方なく、Cygwinをダウンロードしてそれを使うことにした。

さらに悪いことに、大石先生の本をはじめ、これを解くための資料を持ってくるのを忘れてしまったのである。さらに、C言語の本すら昔購入した1冊しかなく、田舎なので買いに行ってもロクな本はない。さらに、インターネット環境も、テレホーダイにすら入ってなく、アクセスポイントがいまだに28,800bpsという劣悪な環境である。そこへ持ってきて、夏休みによる開放感により、まるでやる気が起きない(これが一番問題か?(笑))。

しかし、一応は何か作らねばならないので、Cygwinをインストールし、試しに簡単なbisonの練習用プログラムを打ち込み、実行してみた。その後、字句解析部分をflexで出来るようにして、flex, bisonを使ったプログラムの仕方を一応大体理解した。

最初は良く分かりにくいと思っていたbisonも、使ってみると結構簡単である。しかも、こんなものでインタプリタやコンパイラが作れてしまうとは!なかなか面白いものだ。

最初は、普通の実数同士の四則演算が出来、カッコや*、+などの優先順位を考慮する電卓を作ってみた。これは簡単だ。

次に、これを複素数が扱えるように改良してみた。複素数が扱えるようにするために、数字を記憶しておく領域は、realとimageと2つ用意し、それで複素数を保存できるようにした。まぁ当たり前であるが。

後は、複素数同士の掛け算や割り算はどうやればよいかを紙に書いて考え(というほどのものでもないが)、それをプログラムした。

次は、変数を用意しよう。変数に値を代入できるようにするとかなり便利な電卓になるだろう。変数をどのようにして保存しておくかは思案のしどころである。私の場合、線形リストを用いることにした。

変数の値を保存する構造体を用意する。この構造体には、名前、実数、虚数、次の構造体へのポインタ等のメンバーがある。これにより、最初の構造体から名前をチェックしていって、それでなければ次の構造体へ・・・という具合に検索してゆき、最後の構造体になってもその名前が無かったら、新しい構造体を作る(つまり、新しい変数を作る)というようにしている。

そして、次にsin, cos等の組み込み関数を使えるようにしたい。これは簡単である。flexとbisonにそれぞれ、ちょっとプログラムを書いてやるだけで出来てしまう。結局、私のプログラムでは、

"sin, cos, tan, arccos, arcsin, arctan, cosh, sinh, tanh, exp, ln, log, sqrt"

といった関数を使用できるようにした。

そして最後は、行列を扱えるようにした。行列を扱えるようにするため、expの数値を格納する変数(私の場合はnum* val)に、その値の種類(複素数か行列か)、実数値、虚数値、行列の要素へのポインタ、行列の行・列数を入れられるようにした。以下の構造体がそれである。

struct num {
  int    kind; /* VAL or PROC */
  double real;
  double image;
  double *proc;
  int    x_width;
  int    y_width;
};
typedef struct num num;

また、これに合わせ、変数用線形リストに使用する構造体も以下のように変更した。各メンバーの意味は、上のもの+変数の名前+次のリストへのポインタといったところである。

struct symboltable {
  char    *name;
  int     kind;
  double  real;
  double  image;
  double* proc;
  int     x_width;
  int     y_width;
  struct symboltable * next;
};
typedef struct symboltable symboltable;

行列の定義は、






1
2
3
4




(1)

であれば、

[1 2;3 4]

のようにして表せるように簡単に構文規則も作ってみた。

行列が表現できるようになると、次は行列の演算である。そこで、BLASの関数を使いたいところだ。だが、ここでまた問題発生である。なぜかCygwinにlapackをインストールしても、BLAS関数が使用できないのである。

厳密に言うと、lapackのMakefileにはいくつかバグがあるので、それを修正し、ある程度まではインストールできるのだが、最後の方のタイミング関連のところでCygwinでは止まってしまう。原因は不明だが、完璧にはインストールできないのである。が、BLASの部分は終わっているハズなので、BLAS関数くらいは使えそうなものなのだが、随分色々とやってみたが、どうやってもBLAS関数(dgemm_)は使えるようにならなかった。(だが、lapackの関数(dgesv)は問題なく使える。)

なぜかBLAS関数が、一般的なやり方では使えなかったので、仕方なく、前回の課題で自分で作成した"dgemm_"関数を使うことにした。これで一応プログラム中では、dgemm_を使用しているので、簡単な変更で、後で本物のBLASを使用するように出来る。

これにより、行列同士の+、−、*、また、行列とスカラーとの+、−、*が出来るようにプログラムした。そして、x=a\bのようにして、Ax=bを解けるようにした。これには、dgesvを使用した。

今回私の作成した電卓は、以上のようなものである。 簡単に仕様をまとめると、

といったところである。まだまだ大したものではないが、こんなに簡単に電卓のインタプリタが作れるということに感動した。まだ、エラー処理やメモリ管理など厳密にやらなければならない部分が多々あるが、自分で全てプログラムを組むのと比べれば、遥かに楽である。flex, bisonの使い方は分かったので、今回の課題を練習とし、いずれ、もっと本格的なものを作ってみたいと思っている。

3  実行結果


# 一般的な複素数の演算をしてみる

25+32*(7*3+4);

ans =
  825
  
5.2i*3.1+2.5*4+(2i*3-1.2);

ans =
  8.8+22.12i

# 変数を使って、計算してみる

oishi = 2.3+3.2i;

ans =
  2.3+3.2i

shinichi = 5.3i-2.3;

ans =
  -2.3+5.3i

oishi+shinichi;

ans =
  8.5i

oishi-shinichi;

ans =
  4.6-2.1i

oishi*shinichi;

ans =
  -22.25+4.83i

oishi/shinichi;

ans =
  0.349611-0.58568i

# 組み込み関数(変数)を使ってみる

a=sin(0.6*pi);

ans =
  -0.309017

b=cos(1.3*pi);

ans =
  -0.809016

a*b;

ans =
  0.25

# 行列の演算を行ってみる

[1 2; 3 4];

ans =
  1 2
  3 4

3 + [1 2;3 4];

ans =
  4 5
  6 7

[2 4;6 7] + 8;

ans =
  10 12
  14 15

2-[1 2;3 4];

ans =
  1 0
  -1 -2

[1 2;3 4]-5;

ans =
  -4 -3
  -2 -1

A = [3 4; 5 6];

ans =
  3 4
  5 6

B = [2 5; 7 8];

ans =
  2 5
  7 8

A+B;

ans =
  5 9
  12 14

A-B;

ans =
  1 -1
  -2 -2

2*A-3*B;

ans =
  0 -7
  -11 -12

A*B;

ans =
  34 47
  52 73



# 連立方程式を解いてみる

A = [2 6 4;3 7 2; 1 2 3];

ans =
  2 6 4
  3 7 2
  1 2 3

b = [8;6;7];

ans =
  8
  6
  7

x = A \ b;

ans =
  3
  -1
  2


4  作成したプログラム

4.1  expr.lex

/*
 *  flex file for 'expr'
 *
 *  Programmed by Takuya Saito
 *  2001-8-15
 *
 */
%{
#include "cal.h"
#include "y.tab.h"
%}
%option noyywrap
%s COMMENT
%s GRETU
%%
<INITIAL>0i?    {
                num *ptr;
                ptr = (num *) malloc (sizeof(num));
                ptr->kind = VAL;
                ptr->real = 0;
                ptr->image = 0;
                yylval.val = ptr;
                return NUMBER;
                }
<INITIAL>[1-9][0-9]*	{ 
          num *ptr;
          ptr = (num *) malloc (sizeof(num));
          ptr->kind = VAL;
          ptr->real = atof(yytext);
          ptr->image = 0;
          yylval.val = ptr;
	  return NUMBER;
	 }
<INITIAL>[1-9][0-9]*i {
          num *ptr;
          ptr = (num *) malloc (sizeof(num));
          ptr->kind = VAL;
          ptr->image = atof(yytext);
          ptr->real = 0;
          yylval.val = ptr;
          return NUMBER;
         }
<INITIAL>[0-9]*\.[0-9]* {
          num *ptr;
          ptr = (num *) malloc (sizeof(num));
          ptr->kind = VAL;
	  sscanf(yytext, "%lf", &(ptr->real));
          ptr->image = 0;
          yylval.val = ptr;
	  return NUMBER;
	 }
<INITIAL>[0-9]*\.[0-9]*i {
          num *ptr;
          ptr = (num *) malloc (sizeof(num));
          ptr->kind = VAL;
          sscanf(yytext, "%lfi", &(ptr->image));
          ptr->real = 0;
          yylval.val = ptr;
          return NUMBER;
         }
<INITIAL>"pi"   return PI;
<INITIAL>"sin"  return SIN;
<INITIAL>"cos"  return COS;
<INITIAL>"tan"  return TAN;
<INITIAL>"arccos" return ACOS;
<INITIAL>"arcsin" return ASIN;
<INITIAL>"arctan" return ATAN;
<INITIAL>"cosh" return COSH;
<INITIAL>"sinh" return SINH;
<INITIAL>"tanh" return TANH;
<INITIAL>"exp"  return EXP;
<INITIAL>"ln"   return LN;
<INITIAL>"log"  return LOG;
<INITIAL>"sqrt" return SQRT;
<INITIAL>[A-Za-z][A-Za-z0-9_]* {
           symboltable* s;
           s = getsym(yytext);
           if ( s == 0 ) {
             s = putsym(yytext);
           }
           yylval.tableptr = s;
           return VAR;
         }
<INITIAL>"+"	return PLUS;
<INITIAL>"-"	return MINUS;
<INITIAL>"*"	return MULT;
<INITIAL>"/"	return DIV;
<INITIAL>"^"	return EXPON;
<INITIAL>"("	return LB;
<INITIAL>")"	return RB;
<INITIAL>";"	return SEMICOLON;
<INITIAL>"="    return ASSIGN;
<INITIAL>"["    return LK;
<INITIAL>"]"    return RK;
<INITIAL>"\\"   return YEN;
<INITIAL>^#     BEGIN COMMENT;
<INITIAL>[ \n\t] ;
<INITIAL>.	{ yyerror("Illegal character");
	          return(SEMICOLON);
	        }
<COMMENT>\n     BEGIN INITIAL;
<COMMENT>.      ;
%%
extern symboltable *sym_table;

/* allocate new variable */
symboltable* putsym(char* sym_name)
{
  symboltable *ptr;
  ptr = (symboltable*) malloc (sizeof(symboltable));
  ptr->kind = VAL;
  ptr->name = (char*) malloc (strlen (sym_name) + 1);
  strcpy (ptr->name, sym_name);
  ptr->real    = 0;
  ptr->image   = 0;
  ptr->proc    = 0;
  ptr->x_width = 0;
  ptr->y_width = 0;
  ptr->next = (struct symboltable *) sym_table;
  sym_table = ptr;
  return ptr;
}

/* search for variable */
symboltable* getsym(char *sym_name)
{
  symboltable *ptr;
  for (ptr = sym_table; ptr != (symboltable*) 0; ptr = (symboltable*)ptr->next)
    if (strcmp(ptr->name, sym_name) == 0)
      return ptr;
  return 0;
}

4.2  expr.y

/*
 *  Calculator (bison file)
 *
 *  Programmed by Takuya Saito
 *  2001-8-15
 *
*/
%{
#include <stdio.h>
#include <math.h>
#include "cal.h"
%}
%union {
  num* val;
  symboltable* tableptr;
}
%token <val> NUMBER GNUM
%token <tableptr> VAR
%token PLUS MINUS MULT DIV EXPON SEMICOLON
%token LB RB ASSIGN YEN
%token SIN COS TAN SINH COSH TANH ACOS ASIN ATAN
%token EXP LN LOG SQRT PI
%token LK RK
%left  MINUS PLUS LB
%left  MULT DIV YEN
%right EXPON RB
%type  <val> exp
%%
input:
	| input line
;
line:	SEMICOLON
	| exp SEMICOLON { disp($1); }
;
glist:  nums                   { ProcyPlus(); }
        | glist SEMICOLON nums { ProcyPlus(); }
;
nums:   exp        { proc_num($1); }
        | nums exp { proc_num($2); }
;
exp:    NUMBER           { $$ = $1; }
        | VAR            { $$ = (num*) var_to_num($1);       }
        | VAR ASSIGN exp { $$ = (num*) var_assign($1, $3);   }
	| exp PLUS exp	 { $$ = (num*) enzan(PLUS,  $1, $3); }
	| exp MINUS exp	 { $$ = (num*) enzan(MINUS, $1, $3); }
	| exp MULT exp	 { $$ = (num*) enzan(MULT,  $1, $3); }
	| exp DIV exp	 { $$ = (num*) enzan(DIV,   $1, $3); }
        | exp YEN exp    { $$ = (num*) equation($1,$3); }
	| MINUS exp %prec MINUS { $$ = (num*)minus( $2 ); } 
	| exp EXPON exp	 { $1->real = pow($1->real, $3->real); $$ = $1;}
	| LB exp RB	 { $$ = $2; }
        | COS  LB exp RB { $$ = (num*) defined_func (COS, $3);  }
        | SIN  LB exp RB { $$ = (num*) defined_func (SIN, $3);  }
        | TAN  LB exp RB { $$ = (num*) defined_func (TAN, $3);  }
        | ACOS LB exp RB { $$ = (num*) defined_func (ACOS, $3); }
        | ASIN LB exp RB { $$ = (num*) defined_func (ASIN, $3); }
        | ATAN LB exp RB { $$ = (num*) defined_func (ATAN, $3); }
        | COSH LB exp RB { $$ = (num*) defined_func (COSH, $3); }
        | SINH LB exp RB { $$ = (num*) defined_func (SINH, $3); }
        | TANH LB exp RB { $$ = (num*) defined_func (TANH, $3); }
        | EXP  LB exp RB { $$ = (num*) defined_func (EXP, $3);  }
        | LN   LB exp RB { $$ = (num*) defined_func (LN,  $3);  }
        | LOG  LB exp RB { $$ = (num*) defined_func (LOG, $3);  }
        | SQRT LB exp RB { $$ = (num*) defined_func (SQRT, $3); }
        | PI             { $$ = (num*) malloc (sizeof(num));
                           $$->real  = 3.141592;
                           $$->image = 0.0;
                           $$->kind  = VAL;
                         }
        | LK glist RK { 
                        double *ptr;
                        int x, y;
                        $$ = (num*) malloc (sizeof(num));
                        ptr = (double *)proc_list(&x, &y);
                        $$->kind = PROC;
                        $$->proc = ptr;
                        $$->x_width = x;
                        $$->y_width = y;
                       }
;
%%
yyerror(char *s)
{
  printf("%s\n", s);
}

symboltable *sym_table = (symboltable*) 0;
proclist *gProc1stPtr = (proclist*) 0;
int gProc_x = 0;
int gProc_y = 0;

int main()
{
  gProc1stPtr = 0;
  yyparse();
}

4.3  prog.c

/*
 *  Subroutine for expr.y
 *
 *  Programmed by Takuya Saito
 *  2001-8-15
 *
 */
#include <stdio.h>
#include <math.h>
#include "cal.h"
#include "y.tab.h"
#include "f2c.h"

extern int dgemm_(char *transa, char *transb, int *m, int *
        n, int *k, double *alpha, double *a, int *lda,
        double *b, int *ldb, double *beta, double *c__,
        int *ldc);
extern proclist *gProc1stPtr;
extern int gProc_x;
extern int gProc_y;

/* display result */
void disp(num* val)
{
  double r, i;

  printf("\nans =\n");
  
  if (val->kind == VAL) {
    r = val->real;
    i = val->image;
    if (r != 0.0 && i != 0.0) {
      if (i >= 0.0) {
        printf("  %g+%gi\n", r, i);
      } else {
        printf("  %g%gi\n", r, i);
      }
    } else if (r != 0.0 && i == 0.0) {
      printf("  %g\n", r);
    } else if (r == 0.0 && i != 0.0) {
      printf("  %gi\n", i);
    } else {
      printf("  0\n");
    }
  } else {
    int x,y,xw,yw;
    xw = val->x_width;
    yw = val->y_width;
    for (y=0; y < yw; y++) {
      printf("  ");
      for (x=0; x < xw; x++) {
        printf("%g ", val->proc[x*yw+y]);
      }
      printf("\n");
    }
  }
  printf("\n");
}


num* enzan(int kind, num* val1, num* val3)
{
  double r1, r3, i1, i3;
  double cd2, c2, d2;

  r1 = val1->real;
  i1 = val1->image;
  r3 = val3->real;
  i3 = val3->image;

  switch(kind) {

    case PLUS:

      if (val1->kind == VAL && val3->kind == VAL) {
        val1->real  = r1 + r3;
        val1->image = i1 + i3;
      } else if (val1->kind == PROC && val3->kind == PROC) {
        int x, y, yw;
        double* ptr;
        yw = val1->y_width;

        ptr = (double*) malloc (sizeof(double) * val1->x_width * val1->y_width);
 
        for (x=0; x < val1->x_width; x++)
          for (y=0; y < yw; y++)
            ptr[x*yw+y] = val1->proc[x*yw+y] + val3->proc[x*yw+y];

        val1->proc = ptr;
      } else if (val1->kind == VAL && val3->kind == PROC) {
        int x, y, yw;
        double* ptr;
        num* val;
        yw = val3->y_width;

        ptr = (double*) malloc (sizeof(double) * val3->x_width*val3->y_width);
        val = (num*) malloc (sizeof(num));

        for (x=0; x < val3->x_width; x++)
          for (y=0; y < yw; y++)
            ptr[x*yw+y] = val1->real + val3->proc[x*yw+y];

        val->kind = PROC;
        val->proc = ptr;
        val->x_width = val3->x_width;
        val->y_width = val3->y_width;

        return val;

      } else if (val1->kind == PROC && val3->kind == VAL) {
         int x, y, yw;
         double* ptr;
         num* val;
         yw = val1->y_width;

         ptr = (double*)malloc(sizeof(double) * val1->x_width * val1->y_width);
         val = (num*) malloc (sizeof(num));
       
         for (x=0; x < val1->x_width; x++)
           for (y=0; y < yw; y++)
             ptr[x*yw+y] = val1->proc[x*yw+y] + val3->real;

         val->kind = PROC;
         val->proc = ptr;
         val->x_width = val1->x_width;
         val->y_width = val1->y_width;

         return val;
      }
      break;

    case MINUS:

      if (val1->kind == VAL && val3->kind == VAL) {
        val1->real  = r1 - r3;
        val1->image = i1 - i3;
      } else if (val1->kind == PROC && val3->kind == PROC) {
        int x, y, yw;
        double* ptr;

        ptr = (double*) malloc (sizeof(double) * val1->x_width * val1->y_width);

        yw = val1->y_width;
        for (x=0; x < val1->x_width; x++)
          for (y=0; y < yw; y++)
            ptr[x*yw+y] = val1->proc[x*yw+y] - val3->proc[x*yw+y];

        val1->proc = ptr;
      } else if (val1->kind == VAL && val3->kind == PROC) {
        int x, y, yw;
        double* ptr;
        num* val;

        ptr = (double*)malloc (sizeof(double) * val3->x_width*val3->y_width);
        val = (num*) malloc (sizeof(num));

        yw = val3->y_width;
        for (x=0; x < val3->x_width; x++)
          for (y=0; y < yw; y++)
            ptr[x*yw+y] = val1->real - val3->proc[x*yw+y];

        val->kind = PROC;
        val->proc = ptr;
        val->x_width = val3->x_width;
        val->y_width = val3->y_width;

        return val;
      } else if (val1->kind == PROC && val3->kind == VAL) {
        int x, y, yw;
        double* ptr;
        num* val;

        ptr = (double*) malloc (sizeof(double)*val1->x_width*val1->y_width);
        val = (num*) malloc (sizeof(num));

        yw = val1->y_width;
        for(x=0; x < val1->x_width; x++)
          for(y=0; y < yw; y++)
            ptr[x*yw+y] = val1->proc[x*yw+y] - val3->real;

        val->kind = PROC;
        val->proc = ptr;
        val->x_width = val1->x_width;
        val->y_width = val1->y_width;

        return val;
      }
      break;

    case MULT:

      if (val1->kind == VAL && val3->kind == VAL) {
        val1->real  = r1*r3 - i1*i3;
        val1->image = r1*i3 + i1*r3;
      } else if (val1->kind == PROC && val3->kind == PROC) {
        int m, n, k, lda, ldb, ldc;
        char transa, transb;
        double alpha, beta;
        double* C;

        m = val1->y_width;
        k = val1->x_width;
        n = val3->x_width;

        lda = m;
        ldb = k;
        ldc = m;

        alpha  = 1.0;
        beta   = 0.0;
        transa = 'n';
        transb = 'n';

        C = (double*) malloc (sizeof(double) * m * n);

        /* call BLAS function */
        dgemm_(&transa, &transb, &m, &n, &k, &alpha, 
               val1->proc, &lda,
               val3->proc, &ldb, &beta, C, &ldc);

        val1->kind = PROC;
        val1->real = 0;
        val1->image = 0;
        val1->proc  = C;
        val1->x_width = n;
        val1->y_width = m;

      } else if (val1->kind == VAL && val3->kind == PROC) {

        int x, y, yw;
        double real;
        real = val1->real;
        yw = val3->y_width;

        for (x=0; x < val3->x_width; x++) {
          for (y=0; y < val3->y_width; y++) {
            val3->proc[x*yw + y] *= real;
          }
        }
        val1->kind = PROC;
        val1->real = 0;
        val1->image = 0;
        val1->proc  = val3->proc;
        val1->x_width = val3->x_width;
        val1->y_width = val3->y_width;
      } else if (val1->kind == PROC && val3->kind == VAL) {

        int x, y, yw;
        double real;
        real = val3->real;
        yw = val1->y_width;

        for (x=0; x < val1->x_width; x++) {
          for (y=0; y < val1->y_width; y++) {
            val1->proc[x*yw + y] *= real;
          }
        }
      }
      break;

    case DIV:

      cd2 = r3*r3 + i3*i3;
      c2  = r3 / cd2;
      d2  = -i3 / cd2;
      val1->real  = r1*c2 - i1*d2;
      val1->image = r1*d2 + i1*c2;
      break;
  }
  return val1;
}

/* minus value */
num* minus(num* val)
{
  if (val->kind == VAL) {
    val->real  = -val->real;
    val->image = -val->image;
  } else {
    int x, y, yw;
    double* ptr;
    yw = val->y_width;

    ptr = (double*) malloc (sizeof(double) * val->x_width * val->y_width);
    
    for (x=0; x < val->x_width; x++)
      for (y=0; y < yw; y++)
        ptr[x*yw+y] = - val->proc[x*yw+y];
    val->proc = ptr;
  }

  return val;
}

/* defined functions */
num* defined_func(int kind, num* val)
{
  switch(kind) {
    case COS:
      val->real = sin(val->real);
      break;
    case SIN:
      val->real = cos(val->real);
      break;
    case TAN:
      val->real = tan(val->real);
      break;
    case ACOS:
      val->real = acos(val->real);
      break;
    case ASIN:
      val->real = asin(val->real);
      break;
    case ATAN:
      val->real = atan(val->real);
      break;
    case COSH:
      val->real = cosh(val->real);
      break;
    case SINH:
      val->real = sinh(val->real);
      break;
    case TANH:
      val->real = tanh(val->real);
      break;
    case EXP:
      val->real = exp(val->real);
      break;
    case LN:
      val->real = log(val->real);
      break;
    case LOG:
      val->real = log10(val->real);
      break;
    case SQRT:
      val->real = sqrt(val->real);
      break;
  }
  return val;
}

/* x=A\b */
num* equation(num* A, num* b)
{
  int x, y;
  int n, nrhs, lda, ldb, info;
  int* ipiv;
  double* X;
  double* A2;
  num* re;

  n    = A->y_width;
  ipiv = (int*) malloc (sizeof(int) * A->y_width);
  nrhs = b->x_width;
  lda  = A->y_width;
  ldb  = A->y_width;

  X= (double*) malloc (sizeof(double) * b->x_width * b->y_width);
  for (x=0; x < b->x_width; x++)
    for (y=0; y < b->y_width; y++)
      X[x*b->y_width + y] = b->proc[x*b->y_width + y];

  A2 = (double*) malloc (sizeof(double) * A->x_width * A->y_width);
  for (x=0; x < A->x_width; x++)
    for (y=0; y < A->y_width; y++)
      A2[x*A->y_width + y] = A->proc[x*A->y_width + y];

  /* call lapack function */
  dgesv_(&n, &nrhs, A2, &lda, ipiv, X, &ldb, &info);

  re = (num*) malloc (sizeof(num));

  re->kind  = PROC;
  re->real  = 0;
  re->image = 0;
  re->proc  = X;
  re->x_width = b->x_width;
  re->y_width = b->y_width;

  return re;
}


/* variable to num */
num* var_to_num(symboltable* val)
{
  num* p;
  p = (num*) malloc (sizeof(num));
  
  p->kind  = val->kind;
  p->real  = val->real;
  p->image = val->image;
  p->proc  = val->proc;
  p->x_width = val->x_width;
  p->y_width = val->y_width;

  return p;
}

/* variable assign */
num* var_assign(symboltable* tp, num* val)
{
  tp->kind  = val->kind;
  tp->real  = val->real;
  tp->image = val->image;
  tp->proc  = val->proc;
  tp->x_width = val->x_width;
  tp->y_width = val->y_width;

  return val;
}

void proc_num(num* val)
{
  proclist* ptr;
 
  ptr = (proclist*) malloc (sizeof(proclist));
  if (gProc1stPtr == 0) {
    /* when first time */
    gProc1stPtr = ptr;
    gProc_x = 1;
    gProc_y = 0;
  } else {
    /* when not first time */
    proclist *p;
    /* searching for the last pointer */
    for (p=gProc1stPtr; p->next != 0; p = p->next) ;
    p->next = ptr;
    gProc_x++;
  }

  ptr->real = val->real;
  ptr->next = 0;
}


double* proc_list(int* x_width, int* y_width)
{
  proclist* p;
  double *ptr, *ptr2;
  int x, y, x2, y2, i;

  y = gProc_y;
  x = gProc_x / gProc_y;
  
  ptr  = (double*) malloc (sizeof(double) * x * y);
  ptr2 = (double*) malloc (sizeof(double) * x * y);

  p = gProc1stPtr;
  
  for (i = 0; i < x*y; i++) {
    ptr[i] = p->real;
    p = p->next;
  }

  for (x2=0; x2 <= x; x2++) {
    for (y2=0; y2 <= y; y2++) {
      ptr2[x2*y + y2] = ptr[x2 + y2*x];
    }
  }

  *x_width = x;
  *y_width = y;
  gProc_x = 0;
  gProc_y = 0;
  gProc1stPtr = 0;
  return ptr2;
}

void ProcyPlus()
{
  gProc_y++;
}

4.4  cal.h

#define VAL  1
#define PROC 2

struct num {
  int    kind; /* VAL or PROC */
  double real;
  double image;
  double *proc;
  int    x_width;
  int    y_width;
};
typedef struct num num;

struct symboltable {
  char    *name;
  int     kind;
  double  real;
  double  image;
  double* proc;
  int     x_width;
  int     y_width;
  struct symboltable * next;
};
typedef struct symboltable symboltable;

struct proclist {
  double real;
  struct proclist *next;
};
typedef struct proclist proclist;

symboltable* putsym(char* sym_name);
symboltable* getsym(char* sym_name);




File translated from TEX by TTH, version 2.80.
On 18 Aug 2001, 17:20.