//example of complex linear system solver.  Study gnlinex.cpp before
//studying this example.  The problem is the same except for complex
//numbers.
//The example is the separable kernel m(j,k)=(j+i)*(j+i),r(j)=j+i
//The exact soln is x(k)=(k+i)/(lambda-1) where 
//lambda=(1/3)(size)^3+(1/2)(size)^2-(5/6)(size) + (size)*(size+1) i
//i = sqr(-1) in above notes.

#include <iostream.h>
#include <stdio.h>
#include <stdlib.h>
#include <conio.h>
#include <math.h>

const int precision = 8;               
const int zmaxsize=51;
const int cmaxsize=3;

#include "genprec.h"
#include "cgenprec.h"
#include "apps.h"

int size;

complexmp mat[cmaxsizepad][cmaxsizepad],rhs[cmaxsizepad],check[cmaxsizepad];
complexmp csolution[cmaxsizepad];

main()
  {int i,j;
   size=3;
   cout<<'\n';
   cout<<"the test problem is x(m) + (m+i)*(n+i)*x(n) = (m+i)"<<'\n';
   cout<<"for a 3x3 system"<<'\n'; 
   cout<<'\n';
   cout<<"the matrix factor of x(n) and the right hand side for various m,n are"<<'\n';
   cout<<'\n';
   for (i=1;i<size+1;i++)
     {rhs[i]=cinttomp(i,1);
      cout<<i<<'\n';
      printmp(rhs[i]);
      zmessage="rhs(i)";
      zbreak();
      for (j=1;j<size+1;j++)
        {mat[i][j]=cinttomp(j,1)*rhs[i];
         if (i==j)
           {mat[i][j]=cinttomp(1,0)+mat[i][j];}
         cout<<i<<'\t'<<j<<'\n';
         printmp(mat[i][j]);
         zmessage="mat(i,j)";
         zbreak();
        }
     }

   csolution<csolvlin(size,mat,rhs);
   for (i=1;i<size+1;i++)
     {printmp(csolution[i]);
      zmessage="solution";
      zbreak();
     }

   cout<<"check obtained by multiplying left hand side by solution";
   cout<<'\n';
   for (i=1;i<size+1;i++)
     {check[i]=cinttomp(0,0);     //oops
      for (j=1;j<size+1;j++)
        {check[i]=check[i]+mat[i][j]*csolution[j];}
      printmp(check[i]);
      cout<<'\n';
     }
   return 0;
  }