mplapack-devel — MPACK developer's mailing list

Re: [Mplapack-devel] Can I compute eigenvalues of ill-conditioned matrix less than 10^-156 ?

[Gang Y. <ee...@gm...>]

Hi Maho,

On Fri, Sep 23, 2011 at 4:59 AM, Maho NAKATA <ch...@ma...> wrote:

> Hi Gang,
>
> From: Gang Yan <ee...@gm...>
> Subject: Re: [Mplapack-devel] Can I compute eigenvalues of ill-conditioned
> matrix less than 10^-156 ?
> Date: Fri, 23 Sep 2011 00:55:35 +0800
>
> >   but the smallest eigenvalues as pointed out ~~~~~~~~~ are negative. So
> it
> > is wrong. I think the program can only compute the eigenvalues which are
> > larger than 1e-156 accurately. Do you think so? Or I am wrong for some
> > reasons? Thanks.
>
> in your case, you use
>
> > //initialization of MPFR
> >    int default_prec = 4096;
> >    mpfr_set_default_prec(default_prec);
>
>
 Yes.


> then, machine epsilon is approximately 7.458340731e-155 (*).
> In this case, relative errors of
> eigenvalues can be less accurate than 7.458340731e-155.
> you are doing something complicated in your program and
> I guess some rouding errors included. In my opinion,
> therefore, you found some small negative eigenvalues.
>
> (*) you can check this by inserting following code and rerun your
> program.
>    mpreal a = Rlamch_mpfr("E");
>    mpfr_out_str(stdout, 10, 10, a, MPFR_RNDN);
>    printf("\n");
> Thanks,
>

 I will check it. Many thanks.


> -- Nakata Maho http://accc.riken.jp/maho/ , JA OOO
> http://ja.openoffice.org/
> http://blog.goo.ne.jp/nakatamaho/ ,GPG:
> http://accc.riken.jp/maho/maho.pgp.txt
>

best,
gang

Re: [Mplapack-devel] Can I compute eigenvalues of ill-conditioned matrix less than 10^-156 ?

[Maho N. <ch...@ma...>]

Hi Gang,

From: Gang Yan <ee...@gm...>
Subject: Re: [Mplapack-devel] Can I compute eigenvalues of ill-conditioned matrix less than 10^-156 ?
Date: Fri, 23 Sep 2011 00:55:35 +0800

>   but the smallest eigenvalues as pointed out ~~~~~~~~~ are negative. So it
> is wrong. I think the program can only compute the eigenvalues which are
> larger than 1e-156 accurately. Do you think so? Or I am wrong for some
> reasons? Thanks.

in your case, you use 

> //initialization of MPFR
>    int default_prec = 4096;
>    mpfr_set_default_prec(default_prec);

then, machine epsilon is approximately 7.458340731e-155 (*).
In this case, relative errors of 
eigenvalues can be less accurate than 7.458340731e-155.
you are doing something complicated in your program and
I guess some rouding errors included. In my opinion,
therefore, you found some small negative eigenvalues.

(*) you can check this by inserting following code and rerun your
program.
    mpreal a = Rlamch_mpfr("E");
    mpfr_out_str(stdout, 10, 10, a, MPFR_RNDN);
    printf("\n");
Thanks,
-- Nakata Maho http://accc.riken.jp/maho/ , JA OOO http://ja.openoffice.org/
http://blog.goo.ne.jp/nakatamaho/ ,GPG: http://accc.riken.jp/maho/maho.pgp.txt

Re: [Mplapack-devel] Can I compute eigenvalues of ill-conditioned matrix less than 10^-156 ?

[Gang Y. <ee...@gm...>]

Hi Maho,

   Thanks so much.

On Thu, Sep 22, 2011 at 10:41 PM, Maho NAKATA <ch...@ma...> wrote:

> Hi Gang Yan,
>
> sorry for long delay.
>
> From: Gang Yan <ee...@gm...>
> Subject: Re: [Mplapack-devel] Can I compute eigenvalues of ill-conditioned
> matrix less than 10^-156 ?
> Date: Tue, 13 Sep 2011 01:35:22 -0400
>
> > Many thanks, Maho. Attached please find the source code.
> >
> > A20.dat is the input matrix, eigenvalue_mpfr.cpp is the source code
> modified
> > from your example code.
> ok,
>
> > Place the two files into the folder mpack-0.6.7/examples/mlapack/ and use
> > commands "make" and  "./eigenvalue_mpfr",  it will show the results of
> the
> > largest and smallest eigenvalues.
> done.
>
> > It is shown that when the variable tt (a time variable in my program) is
> > very small, the matrix will become extremely ill-conditioned, and the
> > smallest eigenvalue will be less than 1e-156 and thus become negative.
> following output is correct wrong result?
>

  The output is the same as mine. But it is wrong. Because the matrix is
positive definite and all the eigenvalues should be positive, but ...

> ---
> OK
> [ [ 1.644791426e1]; [ 1.722920009e1]; [ 1.747945946e1]; [ 1.766113939e1]; [
> 1.776131407e1]; [ 1.810655068e1]; [ 1.869304632e1]; [ 1.913033628e1]; [
> 1.941946805e1]; [ 1.992760133e1]; [ 2.000190806e1]; [ 2.041030420e1]; [
> 2.073058106e1]; [ 2.092543236e1]; [ 2.108105957e1]; [ 2.166945660e1]; [
> 2.212463197e1]; [ 2.255747777e1]; [ 2.281199141e1]; [ 2.583112709e1] ]
> OK
> 1.000000000e-4 9.980027172e-5 9.980027105e-5 -3.323972836e-158
>

         ~~~~~~~~~~~~~~~

> 2.000000000e-4 1.992021715e-4 1.992021662e-4 -8.891981975e-159
>

          ~~~~~~~~~~~~~~~

> 4.000000000e-4 3.968173357e-4 3.968172935e-4 -1.165480641e-159
>

         ~~~~~~~~~~~~~~~~

> 8.000000000e-4 7.873381119e-4 7.873377765e-4 3.018299932e-160
> 1.600000000e-3 1.549895802e-3 1.549893166e-3 -2.121601048e-158
>

         ~~~~~~~~~~~~~~~~

  but the smallest eigenvalues as pointed out ~~~~~~~~~ are negative. So it
is wrong. I think the program can only compute the eigenvalues which are
larger than 1e-156 accurately. Do you think so? Or I am wrong for some
reasons? Thanks.
  P.S., the mid two columns are the largest eigenvalues, they all are
correct.


3.200000000e-3 3.003823837e-3 3.003803478e-3 5.848200770e-153
> 6.400000000e-3 5.647598950e-3 5.647447187e-3 3.014126150e-141
> 1.280000000e-2 1.002517109e-2 1.002411737e-2 1.456356388e-129
> 2.560000000e-2 1.606256118e-2 1.605621952e-2 6.184498678e-118
> 5.120000000e-2 2.193695924e-2 2.190825888e-2 2.028549225e-106
> 1.024000000e-1 2.495037779e-2 2.487359809e-2 3.969087087e-95
> 2.048000000e-1 2.545631687e-2 2.535104125e-2 2.761963837e-84
> 4.096000000e-1 2.546875365e-2 2.536157283e-2 2.427716249e-74
> 8.192000000e-1 2.546876239e-2 2.536157891e-2 3.411560868e-66
> 1.638400000 2.546876239e-2 2.536157891e-2 1.315185105e-61
> 3.276800000 2.546876239e-2 2.536157891e-2 2.514747138e-61
> 6.553600000 2.546876239e-2 2.536157891e-2 2.514747138e-61
> 1.310720000e1 2.546876239e-2 2.536157891e-2 2.514747138e-61
> 2.621440000e1 2.546876239e-2 2.536157891e-2 2.514747138e-61
> 5.242880000e1 2.546876239e-2 2.536157891e-2 2.514747138e-61
> ---
>
> thanks
> -- Nakata Maho http://accc.riken.jp/maho/ , JA OOO
> http://ja.openoffice.org/
> http://blog.goo.ne.jp/nakatamaho/ ,GPG:
> http://accc.riken.jp/maho/maho.pgp.txt
>

best,
gang

Re: [Mplapack-devel] Can I compute eigenvalues of ill-conditioned matrix less than 10^-156 ?

[Maho N. <ch...@ma...>]

Hi Gang Yan,

sorry for long delay.

From: Gang Yan <ee...@gm...>
Subject: Re: [Mplapack-devel] Can I compute eigenvalues of ill-conditioned matrix less than 10^-156 ?
Date: Tue, 13 Sep 2011 01:35:22 -0400

> Many thanks, Maho. Attached please find the source code.
> 
> A20.dat is the input matrix, eigenvalue_mpfr.cpp is the source code modified
> from your example code.
ok,

> Place the two files into the folder mpack-0.6.7/examples/mlapack/ and use
> commands "make" and  "./eigenvalue_mpfr",  it will show the results of the
> largest and smallest eigenvalues.
done.

> It is shown that when the variable tt (a time variable in my program) is
> very small, the matrix will become extremely ill-conditioned, and the
> smallest eigenvalue will be less than 1e-156 and thus become negative.
following output is correct wrong result?
---
OK
[ [ 1.644791426e1]; [ 1.722920009e1]; [ 1.747945946e1]; [ 1.766113939e1]; [ 1.776131407e1]; [ 1.810655068e1]; [ 1.869304632e1]; [ 1.913033628e1]; [ 1.941946805e1]; [ 1.992760133e1]; [ 2.000190806e1]; [ 2.041030420e1]; [ 2.073058106e1]; [ 2.092543236e1]; [ 2.108105957e1]; [ 2.166945660e1]; [ 2.212463197e1]; [ 2.255747777e1]; [ 2.281199141e1]; [ 2.583112709e1] ]
OK
1.000000000e-4 9.980027172e-5 9.980027105e-5 -3.323972836e-158 
2.000000000e-4 1.992021715e-4 1.992021662e-4 -8.891981975e-159 
4.000000000e-4 3.968173357e-4 3.968172935e-4 -1.165480641e-159 
8.000000000e-4 7.873381119e-4 7.873377765e-4 3.018299932e-160 
1.600000000e-3 1.549895802e-3 1.549893166e-3 -2.121601048e-158 
3.200000000e-3 3.003823837e-3 3.003803478e-3 5.848200770e-153 
6.400000000e-3 5.647598950e-3 5.647447187e-3 3.014126150e-141 
1.280000000e-2 1.002517109e-2 1.002411737e-2 1.456356388e-129 
2.560000000e-2 1.606256118e-2 1.605621952e-2 6.184498678e-118 
5.120000000e-2 2.193695924e-2 2.190825888e-2 2.028549225e-106 
1.024000000e-1 2.495037779e-2 2.487359809e-2 3.969087087e-95 
2.048000000e-1 2.545631687e-2 2.535104125e-2 2.761963837e-84 
4.096000000e-1 2.546875365e-2 2.536157283e-2 2.427716249e-74 
8.192000000e-1 2.546876239e-2 2.536157891e-2 3.411560868e-66 
1.638400000 2.546876239e-2 2.536157891e-2 1.315185105e-61 
3.276800000 2.546876239e-2 2.536157891e-2 2.514747138e-61 
6.553600000 2.546876239e-2 2.536157891e-2 2.514747138e-61 
1.310720000e1 2.546876239e-2 2.536157891e-2 2.514747138e-61 
2.621440000e1 2.546876239e-2 2.536157891e-2 2.514747138e-61 
5.242880000e1 2.546876239e-2 2.536157891e-2 2.514747138e-61 
---

thanks
-- Nakata Maho http://accc.riken.jp/maho/ , JA OOO http://ja.openoffice.org/
http://blog.goo.ne.jp/nakatamaho/ ,GPG: http://accc.riken.jp/maho/maho.pgp.txt

Re: [Mplapack-devel] Can I compute eigenvalues of ill-conditioned matrix less than 10^-156 ?

[Gang Y. <ee...@gm...>]

Many thanks, Maho. Attached please find the source code.

A20.dat is the input matrix, eigenvalue_mpfr.cpp is the source code modified
from your example code.

Place the two files into the folder mpack-0.6.7/examples/mlapack/ and use
commands "make" and  "./eigenvalue_mpfr",  it will show the results of the
largest and smallest eigenvalues.
It is shown that when the variable tt (a time variable in my program) is
very small, the matrix will become extremely ill-conditioned, and the
smallest eigenvalue will be less than 1e-156 and thus become negative.

thanks,
gang

On Mon, Sep 12, 2011 at 11:53 PM, Maho NAKATA <ch...@ma...> wrote:

> From: Gang Yan <ee...@gm...>
> Subject: Re: [Mplapack-devel] Can I compute eigenvalues of ill-conditioned
> matrix less than 10^-156 ?
> Date: Tue, 13 Sep 2011 11:45:47 +0800
>
> > Hi, Maho,
> >
> >
> > Gang, if you want to use MPACK, as Kouya-sensei wrote, please use
> >> GMP or MPFR version (to enlarge your precision).
> >>
> > I published a tutorial in Japanese (sorry), and
> >> I'll upload to the web page when I can publish my original one at the
> web.
> >
> >
> >> In the mean time, if you can use GMP version and setting
> >> MPACK_GMP_PRECISION environment variable to change pricision
> dynamically.
> >>
> >
> >
> > Yes, I used the GMP and MPFR version and tuned the precision variable
> from
> > 256 to 2048, but the results are truncated at ~10^-156. When the
> eigenvalues
> > are smaller than ~10^-156, they will become negative. How can I improve
> the
> > precision beyond 10^-156? Thanks.
>
> Hm, could you please check your inputs, so that not to include trancated
> values.
> or could you please send me your source codes.
>
> thanks
>  Nakata Maho
>

Re: [Mplapack-devel] Can I compute eigenvalues of ill-conditioned matrix less than 10^-156 ?

[Maho N. <ch...@ma...>]

From: Gang Yan <ee...@gm...>
Subject: Re: [Mplapack-devel] Can I compute eigenvalues of ill-conditioned matrix less than 10^-156 ?
Date: Tue, 13 Sep 2011 11:45:47 +0800

> Hi, Maho,
> 
> 
> Gang, if you want to use MPACK, as Kouya-sensei wrote, please use
>> GMP or MPFR version (to enlarge your precision).
>>
> I published a tutorial in Japanese (sorry), and
>> I'll upload to the web page when I can publish my original one at the web.
> 
> 
>> In the mean time, if you can use GMP version and setting
>> MPACK_GMP_PRECISION environment variable to change pricision dynamically.
>>
> 
> 
> Yes, I used the GMP and MPFR version and tuned the precision variable from
> 256 to 2048, but the results are truncated at ~10^-156. When the eigenvalues
> are smaller than ~10^-156, they will become negative. How can I improve the
> precision beyond 10^-156? Thanks.

Hm, could you please check your inputs, so that not to include trancated
values.
or could you please send me your source codes.

thanks
 Nakata Maho

Re: [Mplapack-devel] Can I compute eigenvalues of ill-conditioned matrix less than 10^-156 ?

[Gang Y. <ee...@gm...>]

Hi, Maho,


Gang, if you want to use MPACK, as Kouya-sensei wrote, please use
> GMP or MPFR version (to enlarge your precision).
>
I published a tutorial in Japanese (sorry), and
> I'll upload to the web page when I can publish my original one at the web.


> In the mean time, if you can use GMP version and setting
> MPACK_GMP_PRECISION environment variable to change pricision dynamically.
>


Yes, I used the GMP and MPFR version and tuned the precision variable from
256 to 2048, but the results are truncated at ~10^-156. When the eigenvalues
are smaller than ~10^-156, they will become negative. How can I improve the
precision beyond 10^-156? Thanks.


>
> Thanks
>  Nakata Maho
>

best,
gang

Re: [Mplapack-devel] Can I compute eigenvalues of ill-conditioned matrix less than 10^-156 ?

[Gang Y. <ee...@gm...>]

Hi,  Tomonori,

 Thank you very much. Yes, for very ill-conditioned matrices, eigenvalues
are very sensitive to round-off errors when computing with float-point. I
think I can obtain the smallest eigenvalues if using higher precision
arithmetic.

best,
gang

On Mon, Sep 12, 2011 at 11:01 PM, Tomonori Kouya <ju...@qu...>wrote:

> Hi, Gang. I have tried to solve some eigenvalue problems using MPFR/GMP.
>
> Very ill-conditionded matrices such as Hilbert or Lotkin matrices, are
> too sensitive to keep their theoretical properties due to initial error
> or round-off error in computing process. In case of Hilbert matrix, it
> is theoretically positive definite but some of smallest eigenvalues can
> be minus due to initial error in approximated elements of Hilbert matrix.
>
> I recommend to know how sensitive your matrix is for such error by using
> variable precision of MPFR/GMP in order to solve your troubles.
>
> (2011/09/12 23:35), Gang Yan wrote:
> >     Now I have a problem when I calculate the eigenvalues of an extremely
> > ill-conditioned matrix, in a problem of control system research. The
> matrix
> > is not very large, e.g., 30X30, but extremely ill-conditioned, i.e.,
> having
> > the largest eigenvalue ~10^-2 and the smallest eigenvalue>  10^-156. I
> found
> > that the eigenvalue programme using GMP or MPFR can only compute
> eigenvalues
> > <10^-156, or I will get some minus eigenvalues which is not wrong because
> > the matrix is positive-definition. I use a 64-bit notebook.
>
>
> ------------------------------------------------------------------------------
> Doing More with Less: The Next Generation Virtual Desktop
> What are the key obstacles that have prevented many mid-market businesses
> from deploying virtual desktops?   How do next-generation virtual desktops
> provide companies an easier-to-deploy, easier-to-manage and more affordable
> virtual desktop model.http://www.accelacomm.com/jaw/sfnl/114/51426474/
> _______________________________________________
> Mplapack-devel mailing list
> Mpl...@li...
> https://lists.sourceforge.net/lists/listinfo/mplapack-devel
>

Re: [Mplapack-devel] Can I compute eigenvalues of ill-conditioned matrix less than 10^-156 ?

[Maho N. <ch...@ma...>]

Hi Gang and Kouya-sensei

Gang, if you want to use MPACK, as Kouya-sensei wrote, please use
GMP or MPFR version (to enlarge your precision).
I published a tutorial in Japanese (sorry), and
I'll upload to the web page when I can publish my original one at the web.

In the mean time, if you can use GMP version and setting
MPACK_GMP_PRECISION environment variable to change pricision dynamically.

Thanks
 Nakata Maho



From: Tomonori Kouya <ju...@qu...>
Subject: Re: [Mplapack-devel] Can I compute eigenvalues of ill-conditioned matrix less than 10^-156 ?
Date: Tue, 13 Sep 2011 00:01:53 +0900

> Hi, Gang. I have tried to solve some eigenvalue problems using MPFR/GMP.
> 
> Very ill-conditionded matrices such as Hilbert or Lotkin matrices, are 
> too sensitive to keep their theoretical properties due to initial error 
> or round-off error in computing process. In case of Hilbert matrix, it 
> is theoretically positive definite but some of smallest eigenvalues can 
> be minus due to initial error in approximated elements of Hilbert matrix.
> 
> I recommend to know how sensitive your matrix is for such error by using 
> variable precision of MPFR/GMP in order to solve your troubles.
> 
> (2011/09/12 23:35), Gang Yan wrote:
>>     Now I have a problem when I calculate the eigenvalues of an extremely
>> ill-conditioned matrix, in a problem of control system research. The matrix
>> is not very large, e.g., 30X30, but extremely ill-conditioned, i.e., having
>> the largest eigenvalue ~10^-2 and the smallest eigenvalue>  10^-156. I found
>> that the eigenvalue programme using GMP or MPFR can only compute eigenvalues
>> <10^-156, or I will get some minus eigenvalues which is not wrong because
>> the matrix is positive-definition. I use a 64-bit notebook.
> 
> ------------------------------------------------------------------------------
> Doing More with Less: The Next Generation Virtual Desktop 
> What are the key obstacles that have prevented many mid-market businesses
> from deploying virtual desktops?   How do next-generation virtual desktops
> provide companies an easier-to-deploy, easier-to-manage and more affordable
> virtual desktop model.http://www.accelacomm.com/jaw/sfnl/114/51426474/
> _______________________________________________
> Mplapack-devel mailing list
> Mpl...@li...
> https://lists.sourceforge.net/lists/listinfo/mplapack-devel
> 

Re: [Mplapack-devel] Can I compute eigenvalues of ill-conditioned matrix less than 10^-156 ?

[Tomonori K. <ju...@qu...>]

Hi, Gang. I have tried to solve some eigenvalue problems using MPFR/GMP.

Very ill-conditionded matrices such as Hilbert or Lotkin matrices, are 
too sensitive to keep their theoretical properties due to initial error 
or round-off error in computing process. In case of Hilbert matrix, it 
is theoretically positive definite but some of smallest eigenvalues can 
be minus due to initial error in approximated elements of Hilbert matrix.

I recommend to know how sensitive your matrix is for such error by using 
variable precision of MPFR/GMP in order to solve your troubles.

(2011/09/12 23:35), Gang Yan wrote:
>     Now I have a problem when I calculate the eigenvalues of an extremely
> ill-conditioned matrix, in a problem of control system research. The matrix
> is not very large, e.g., 30X30, but extremely ill-conditioned, i.e., having
> the largest eigenvalue ~10^-2 and the smallest eigenvalue>  10^-156. I found
> that the eigenvalue programme using GMP or MPFR can only compute eigenvalues
> <10^-156, or I will get some minus eigenvalues which is not wrong because
> the matrix is positive-definition. I use a 64-bit notebook.

[Mplapack-devel] Can I compute eigenvalues of ill-conditioned matrix less than 10^-156 ?

[Gang Y. <ee...@gm...>]

Dear Folks,

   Many thanks to Mpack. It is really helpful for computing ill-conditioned
problems.

   Now I have a problem when I calculate the eigenvalues of an extremely
ill-conditioned matrix, in a problem of control system research. The matrix
is not very large, e.g., 30X30, but extremely ill-conditioned, i.e., having
the largest eigenvalue ~10^-2 and the smallest eigenvalue > 10^-156. I found
that the eigenvalue programme using GMP or MPFR can only compute eigenvalues
<10^-156, or I will get some minus eigenvalues which is not wrong because
the matrix is positive-definition. I use a 64-bit notebook.

   Any suggestions or comments? Thanks in advance.

best,
gang

Re: [Mplapack-devel] FORTRAN interface for MPACK

[Maho N. <ch...@ma...>]

Hi Fletcher,

Thanks for your great feedback. I'll provide when I have some time.

From: "Fletcher, John P" <j.p...@as...>
Subject: FORTRAN interface for MPACK
Date: Wed, 06 Apr 2011 16:14:18 +0100

> Nakata Maho
> 
> Here is an example of a FORTRAN interface for MPACK, with some extracts from the makefile needed below.  The C++ file provides interfaces so that the program will link.
> 
> # Extras for Fortran for QD
> FCFLAGS=-I/usr/local/lib/qd
> FCLIBS=-L/usr/local/lib -lqdmod -lqd -lrt -lm -L/usr/lib/gcc/x86_64-linux-gnu/4.4.3 -L/usr/lib/gcc/x86_64-linux-gnu/4.4.3/../../../../lib -L/lib/../lib -L/usr/lib/../lib -L/usr/lib/gcc/x86_64-linux-gnu/4.4.3/../../.. -L/usr/lib/x86_64-linux-gnu -lrt -lgfortran -lm
> FCMAIN=-L/usr/local/lib -lqd_f_main
> 
> %.o: %.f90
>       $(FC) -c $(FCFLAGS) $<
> 
> # QD and MPACK interfaced from FORTRAN
> solve_qd_f: solve_qd_f.o mpack_qd_cpp.o
>       $(CXX) -o $@ solve_qd_f.o mpack_qd_cpp.o \
>         $(LIBFLAGS) $(QDLIBS) $(FCLIBS) $(FCMAIN)
> 
> I  hope this is of interest.
> 
> John Fletcher
> 
> 
> Dr John P. Fletcher Tel: (44) 121 204 3389 (direct line), FAX: (44) 121 204 3678
> Chemical Engineering and Applied Chemistry (CEAC),
> Associate Dean - External Relations,
> School of Engineering and Applied Science (EAS),
> Aston University, Aston Triangle, BIRMINGHAM B4 7ET  U.K.
> 
> 

Re: [Mplapack-devel] Bug in Rgerfs

[Maho N. <ch...@ma...>]

Hi Fletcher,

committed, thanks!

Best regards,
 Nakata Maho

From: "Fletcher, John P" <j.p...@as...>
Subject: Bug in Rgerfs
Date: Wed, 06 Apr 2011 09:18:41 +0100

> Nakata Maho
> 
> I have several things to send you and will start with this one.
> 
> I have been testing Rgerfs for double and dd and did not get the results I get with LAPACK.
> 
> I think the reason is that some of the accesses in Rgerfs to ipiv, B and X start at [1] and not [0], and also access to work is wrongly indexed so that one more variable is accessed at the end.
> 
> I attach two test routines.  solve_double_refine.cpp shows it not working and solve_dd_refine.cpp shows reindexing the ipiv, B and X arrays - I copied them to new arrays - and extra work.  This works except that the error results are wrong - file attached.
> 
> The other file is the correct results for the same problem using the routine from LAPACK.
> 
> I hope this helps.
> 
> John Fletcher
> 
> P.S.  I will send other things later.
> 
> Dr John P. Fletcher Tel: (44) 121 204 3389 (direct line), FAX: (44) 121 204 3678
> Chemical Engineering and Applied Chemistry (CEAC),
> Associate Dean - External Relations,
> School of Engineering and Applied Science (EAS),
> Aston University, Aston Triangle, BIRMINGHAM B4 7ET  U.K.
> 
> 

Re: [Mplapack-devel] Operator<< for MPFR complex

[Maho N. <ch...@ma...>]

Hi Fletcher,

Many thanks for your feedback.
Do you need a commit right of mpack project?
thanks,
Nakata Maho

From: "Fletcher, John P" <j.p...@as...>
Subject: Operator<< for MPFR complex
Date: Wed, 06 Apr 2011 16:05:40 +0100

> Nakata Maho
> 
> The operator<< for mpcomplex is modeled on the one for mpreal.

Opps, that's right :-)

>   std::ostream& operator<<(std::ostream& os, const mpcomplex& v)
> {
>   return os << " ( " << v.real() << " " << v.imag() << " ) ";
> }
> 
> I have attached modified versions of the files mpcomplex.h and mpcomplex.cpp. It provides for output as in the example in trial_mpcomplex.cpp which also documents what I was saying about the libraries.
> 
> Best wishes
> 
> John Fletcher
> 
> Dr John P. Fletcher Tel: (44) 121 204 3389 (direct line), FAX: (44) 121 204 3678
> Chemical Engineering and Applied Chemistry (CEAC),
> Associate Dean - External Relations,
> School of Engineering and Applied Science (EAS),
> Aston University, Aston Triangle, BIRMINGHAM B4 7ET  U.K.
> 
> 

[Mplapack-devel] FORTRAN interface for MPACK

[Fletcher, J. P <j.p...@as...>]

// mpack_qd_cpp.cpp

// Interface qd_complex sqrt function.

// John Fletcher   J.P...@as...  April 2011

#include <mblas_qd.h>
#include <mlapack_qd.h>

// The function bodies to be called from FORTRAN have to be extern "C".
// They also need a trailing _

extern "C" {


  qd_complex qdc_sqrt_(const qd_complex &z)
  {
    return sqrt(z);
  }

  void rgetrf_(const mpackint &m, const mpackint &n, qd_real * A, const mpackint &lda, mpackint *ipiv, mpackint *info )
  {
    Rgetrf(m,n,A,lda,ipiv,info);
  }

  void rgetrs_(const char *trans, const mpackint &n, const mpackint &nrhs, qd_real * A, const mpackint &lda, mpackint *ipiv, qd_real *B, const mpackint &ldb, mpackint *info )
  {
    Rgetrs(trans,n,nrhs,A,lda,ipiv,B,ldb,info);
  }

  void cgetrf_(const mpackint &m, const mpackint &n, qd_complex * A, const mpackint &lda, mpackint *ipiv, mpackint *info )
  {
    Cgetrf(m,n,A,lda,ipiv,info);
  }

  void cgetrs_(const char *trans, const mpackint &n, const mpackint &nrhs, qd_complex * A, const mpackint &lda, mpackint *ipiv, qd_complex *B, const mpackint &ldb, mpackint *info )
  {
    Cgetrs(trans,n,nrhs,A,lda,ipiv,B,ldb,info);
  }

} /* extern "C" */

[Mplapack-devel] Operator<< for MPFR complex

[Fletcher, J. P <j.p...@as...>]

/*
 * Copyright (c) 2010
 *	Nakata, Maho
 * 	All rights reserved.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions
 * are met:
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
 * SUCH DAMAGE.
 *
 */
/*
Complex class declare for the MPFR
*/

#ifndef __MP_COMPLEX_H__
#define __MP_COMPLEX_H__

#include "mpreal.h"
#include <complex>
#include "mpc.h"

#if defined ___MPACK_BUILD_WITH_GMP___
#include "gmpxx.h"
#include "mpc_class.h"
#endif
#if defined ___MPACK_BUILD_WITH_QD___
#include "qd_complex.h"
#endif
#if defined ___MPACK_BUILD_WITH_DD___
#include "dd_complex.h"
#endif

namespace mpfr {
class mpcomplex {
  private:
   mpc_t mpc;

  public:
  static mpc_rnd_t default_rnd;
  static mp_prec_t default_real_prec;
  static mp_prec_t default_imag_prec;
  static int default_base;
  static int double_bits;

//constructor & deconstructor
    mpcomplex();
    mpcomplex(const mpc_t a);
    mpcomplex(const mpfr_t a, const mpfr_t b);
    mpcomplex(const mpf_t a, const mpf_t b);
    mpcomplex(const char *s, mp_prec_t pr = default_real_prec, mp_prec_t pi = default_imag_prec, mpc_rnd_t mode = default_rnd);
    mpcomplex(const mpcomplex& a);
    mpcomplex(const std::complex<double>& a, mp_prec_t pr = default_real_prec, mp_prec_t pi = default_imag_prec, mpc_rnd_t mode = default_rnd);
    mpcomplex(const mpreal& a, const mpreal& b);
    mpcomplex(const double& a, const double& b, mp_prec_t pr = default_real_prec, mp_prec_t pi = default_imag_prec, mpc_rnd_t mode = default_rnd);
    mpcomplex(const char *s, const char *t, mp_prec_t pr = default_real_prec, mp_prec_t pi = default_imag_prec, mpc_rnd_t mode = default_rnd);

    mpcomplex(const mpreal & a);
    mpcomplex(const mpfr_t a);
    mpcomplex(const mpf_t a);
    mpcomplex(const double a, mp_prec_t pr = default_real_prec, mp_prec_t pi = default_imag_prec, mpc_rnd_t mode = default_rnd);
    ~mpcomplex(); 

    mpcomplex& operator=(const mpcomplex& a);
    mpcomplex& operator=(const mpc_t a);
    mpcomplex& operator=(const std::complex<double> a);
    mpcomplex& operator=(const char* s);

    mpcomplex& operator=(const mpreal& a);
    mpcomplex& operator=(const mpfr_t a);
    mpcomplex& operator=(const double a);

    //+
    mpcomplex& operator+=(const mpcomplex& a);
    mpcomplex& operator+=(const mpc_t a);
    mpcomplex& operator+=(const std::complex<double> a);

    mpcomplex& operator+=(const mpreal& a);
    mpcomplex& operator+=(const mpfr_t a);
    mpcomplex& operator+=(const double a);
    const mpcomplex operator+() const;
    mpcomplex& operator++();
    const mpcomplex operator ++ (int);

    //-
    mpcomplex& operator-=(const mpcomplex& a);
    mpcomplex& operator-=(const mpc_t a);
    mpcomplex& operator-=(const std::complex<double> a);

    mpcomplex& operator-=(const mpreal& a);
    mpcomplex& operator-=(const mpfr_t a);
    mpcomplex& operator-=(const double a);
    const mpcomplex operator-() const;
    mpcomplex& operator--();
    const mpcomplex operator -- (int);

    //*
    mpcomplex& operator*=(const mpcomplex& a);
    mpcomplex& operator*=(const mpc_t a);
    mpcomplex& operator*=(const std::complex<double> a);
    mpcomplex& operator*=(const mpreal& a);
    mpcomplex& operator*=(const mpfr_t a);
    mpcomplex& operator*=(const double a);

    // /
    mpcomplex& operator/=(const mpcomplex& a);
    mpcomplex& operator/=(const mpc_t a);
    mpcomplex& operator/=(const std::complex<double> a);
    mpcomplex& operator/=(const mpreal& a);
    mpcomplex& operator/=(const mpfr_t a);
    mpcomplex& operator/=(const double a);

    // Output/ Input
    friend std::ostream& operator<<(std::ostream& os, const mpcomplex& v);
    //    friend std::istream& operator>>(std::istream& is, mpcomplex& v);

    //comparison
    friend bool operator == (const mpcomplex& a, const mpcomplex& b);
    friend bool operator == (const mpcomplex& a, const mpreal& b);
    friend bool operator == (const mpreal& a, const mpcomplex& b);
    friend bool operator != (const mpcomplex& a, const mpcomplex& b);
    friend bool operator != (const mpreal& a, const mpcomplex& b);
    friend bool operator != (const mpcomplex& a, const mpreal& b);

    //random
    friend const mpcomplex urandom_c (gmp_randstate_t& state);

//extraction of real and imaginary parts
    const mpreal real() const
    {
        return mpreal(mpc_realref(mpc));
    }
    const mpreal imag() const
    {
        return mpreal(mpc_imagref(mpc));
    }
    mpreal real()
    {
        return mpreal(mpc_realref(mpc));
    }
    mpreal imag()
    {
        return mpreal(mpc_imagref(mpc));
    }
    //other functions
    friend const mpreal abs(const mpcomplex& a, mpfr_rnd_t mode = mpreal::default_rnd);
    friend const mpreal norm(const mpcomplex& a, mpfr_rnd_t mode = mpreal::default_rnd);
    friend const mpcomplex conj(const mpcomplex& a, mpc_rnd_t mode = default_rnd);
    friend const mpreal arg(const mpcomplex& a, mpfr_rnd_t mode = mpreal::default_rnd);
    friend const mpcomplex proj(const mpcomplex& a, mpc_rnd_t mode = mpreal::default_rnd);

    //Powerfunction and Logarithm    
    friend const mpcomplex sqrt(const mpcomplex& a, mpc_rnd_t mode = default_rnd);
    friend const mpcomplex pow(const mpcomplex& a, const mpcomplex& b,mpc_rnd_t mode = default_rnd);
    friend const mpcomplex exp(const mpcomplex& a, mpc_rnd_t mode = default_rnd);
    friend const mpcomplex log(const mpcomplex& a, mpc_rnd_t mode = default_rnd);

    //trigonometric functions
    friend const mpcomplex sin(const mpcomplex& a, mpc_rnd_t mode = default_rnd);
    friend const mpcomplex cos(const mpcomplex& a, mpc_rnd_t mode = default_rnd);
    friend const mpcomplex tan(const mpcomplex& a, mpc_rnd_t mode = default_rnd);
    friend const mpcomplex sinh(const mpcomplex& a, mpc_rnd_t mode = default_rnd);
    friend const mpcomplex cosh(const mpcomplex& a, mpc_rnd_t mode = default_rnd);
    friend const mpcomplex tanh(const mpcomplex& a, mpc_rnd_t mode = default_rnd);
    friend const mpcomplex asin(const mpcomplex& a, mpc_rnd_t mode = default_rnd);
    friend const mpcomplex acos(const mpcomplex& a, mpc_rnd_t mode = default_rnd);
    friend const mpcomplex atan(const mpcomplex& a, mpc_rnd_t mode = default_rnd);
    friend const mpcomplex asinh(const mpcomplex& a, mpc_rnd_t mode = default_rnd);
    friend const mpcomplex acosh(const mpcomplex& a, mpc_rnd_t mode = default_rnd);
    friend const mpcomplex atanh(const mpcomplex& a, mpc_rnd_t mode = default_rnd);

// Type Conversion operators
    operator std::complex<double>() const;
    operator std::string() const;
    std::string to_string(size_t n = 0, int b = default_base, mpc_rnd_t mode = default_rnd) const;

// Set/Get instance properties
   mp_prec_t get_prec() const;
   mp_prec_t get_prec_re() const;
   mp_prec_t get_prec_im() const;
   void set_prec(mp_prec_t prec, mpc_rnd_t rnd_mode = default_rnd);
   void set_prec2(mp_prec_t pr, mp_prec_t pi, mpc_rnd_t rnd_mode = default_rnd);

#if defined ___MPACK_BUILD_WITH_GMP___
mpcomplex(const mpc_class& a);
operator mpc_class() const;
#endif
#if defined ___MPACK_BUILD_WITH_QD___
mpcomplex(const qd_complex& a);
operator qd_complex() const;
mpcomplex& operator=(const qd_complex& a);
#endif
#if defined ___MPACK_BUILD_WITH_DD___
mpcomplex(const dd_complex& a);
operator dd_complex() const;
mpcomplex& operator=(const dd_complex& a);
#endif
};

//+ addition
const mpcomplex operator+(const mpcomplex& a, const mpcomplex& b);
const mpcomplex operator+(const mpcomplex& a, const std::complex<double> b);
const mpcomplex operator+(const mpcomplex& a, const char* b);

const mpcomplex operator+(const std::complex<double> a, const mpcomplex& b);
const mpcomplex operator+(const char* a, const mpcomplex& b);

const mpcomplex operator+(const mpcomplex& a, const mpreal b);
const mpcomplex operator+(const mpcomplex& a, const double b);
const mpcomplex operator+(const mpcomplex& a, const int b);

const mpcomplex operator+(const mpreal a, const mpcomplex& b);
const mpcomplex operator+(const double a, const mpcomplex& b);
const mpcomplex operator+(const int a, const mpcomplex& b);
//- subtraction
const mpcomplex operator-(const mpcomplex& a, const mpcomplex& b);
const mpcomplex operator-(const mpcomplex& a, const std::complex<double>& b);
const mpcomplex operator-(const mpcomplex& a, const char* b);

const mpcomplex operator-(const std::complex<double>& a, const mpcomplex& b);
const mpcomplex operator-(const char* a, const mpcomplex& b);

const mpcomplex operator-(const mpcomplex& a, const mpreal b);
const mpcomplex operator-(const mpcomplex& a, const double b);
const mpcomplex operator-(const mpcomplex& a, const int b);

const mpcomplex operator-(const mpreal a, const mpcomplex& b);
const mpcomplex operator-(const double a, const mpcomplex& b);
const mpcomplex operator-(const int a, const mpcomplex& b);

//* multiplication
const mpcomplex operator*(const mpcomplex& a, const mpcomplex& b);
const mpcomplex operator*(const mpcomplex& a, const mpreal& b);
const mpcomplex operator*(const mpreal& a, const mpcomplex& b);

/// division
const mpcomplex operator/(const mpcomplex& a, const mpcomplex& b);

inline mp_prec_t mpcomplex::get_prec() const
{
	return mpc_get_prec(mpc);
}

inline mp_prec_t mpcomplex::get_prec_re() const
{
      mp_prec_t pr, pi;
      mpc_get_prec2(&pr, &pi, mpc);
      return pr;
}

inline mp_prec_t mpcomplex::get_prec_im() const
{
      mp_prec_t pr, pi;
      mpc_get_prec2(&pr, &pi, mpc);
      return pi;
}

inline void mpcomplex::set_prec(mp_prec_t prec, mpc_rnd_t rnd_mode)
{
    mpcomplex tmp (*this);
    mpc_init2(mpc, prec);
    mpc_set(mpc, tmp.mpc, default_rnd);
}

inline void mpcomplex::set_prec2(mp_prec_t pr, mp_prec_t pi, mpc_rnd_t rnd_mode)
{
    mpcomplex tmp (*this);
    mpc_init3(mpc, pr, pi);
    mpc_set(mpc, tmp.mpc, default_rnd);
}

//constructor and deconstructor
inline mpcomplex::mpcomplex()
{
   mpc_init3(mpc, default_real_prec, default_imag_prec);
   mpc_set_ui(mpc,0UL,default_rnd);
}

inline mpcomplex::mpcomplex(const mpc_t a) 
{
   mp_prec_t pr, pi;
   mpc_get_prec2(&pr, &pi, a);
   mpc_init3(mpc, pr, pi);
   mpc_set(mpc, a, default_rnd);
}

inline mpcomplex::mpcomplex(const mpfr_t a, const mpfr_t b)
{
   mp_prec_t pr, pi;
   pr = mpfr_get_prec(a);  pi = mpfr_get_prec(b);
   mpc_init3(mpc, pr, pi);
   mpc_set_fr_fr(mpc, a, b, default_rnd);
}

inline mpcomplex::mpcomplex(const mpf_t a, const mpf_t b)
{
   mp_prec_t pr, pi;
   pr = mpf_get_prec(a);  pi = mpf_get_prec(b);
   mpc_init3(mpc, pr, pi);
   mpc_set_f_f(mpc, a, b, default_rnd);
}

inline mpcomplex::mpcomplex(const char *s, mp_prec_t pr, mp_prec_t pi, mpc_rnd_t mode)
{
    mpc_init3(mpc, pr, pi);
    mpc_set_str(mpc, (char *)s, default_base, mode); 
}

inline mpcomplex::mpcomplex(const mpcomplex& a) 
{
   mp_prec_t pr, pi;
   mpc_get_prec2(&pr, &pi, a.mpc);
   mpc_init3(mpc, pr, pi);
   mpc_set(mpc, a.mpc, default_rnd);
}

inline mpcomplex::mpcomplex(const std::complex<double>& a, mp_prec_t pr, mp_prec_t pi, mpc_rnd_t mode)
{
   mpc_init3(mpc, pr, pi);
   mpc_set_d_d(mpc, a.real(), a.imag(), default_rnd);
}

inline mpcomplex::mpcomplex(const mpreal& a, const mpreal& b)
{
   mp_prec_t pr, pi;
   mpreal tmp1(a), tmp2(b);
   pr = a.get_prec();  pi = b.get_prec();
   mpc_init3(mpc, pr, pi);
   mpc_set_fr_fr(mpc, (mpfr_ptr)tmp1, (mpfr_ptr)tmp2, default_rnd);
}

inline mpcomplex::mpcomplex(const double& a, const double& b, mp_prec_t pr, mp_prec_t pi, mpc_rnd_t mode)
{
   mpc_init3(mpc, default_real_prec, default_imag_prec);
   mpc_set_d_d(mpc, a, b,default_rnd);
}

inline mpcomplex::mpcomplex(const char *s, const char *t, mp_prec_t pr, mp_prec_t pi, mpc_rnd_t mode)
{
   mpfr_t a, b;
   mpfr_init2(a, pr); mpfr_init2(b, pi);
   mpfr_set_str(a, s, default_base, MPC_RND_RE(mode));
   mpfr_set_str(b, t, default_base, MPC_RND_IM(mode));
   mpc_init3(mpc, pr, pi);
   mpc_set_fr_fr(mpc, a, b, default_rnd);
   mpfr_clear(a); mpfr_clear(b);
}

inline mpcomplex::mpcomplex(const mpreal& a)
{
   mpreal tmp(a);
   mpc_init3(mpc, default_real_prec, default_imag_prec);
   mpc_set_fr(mpc, (mpfr_ptr)(tmp), default_rnd);
}

inline mpcomplex::mpcomplex(const mpfr_t a)
{
   mp_prec_t pr;
   pr = mpfr_get_prec(a);
   mpc_init2(mpc, pr);
   mpc_set_fr(mpc, a, default_rnd);
}

inline mpcomplex::mpcomplex(const mpf_t a)
{
   mp_prec_t pr;
   pr = mpf_get_prec(a);
   mpc_init2(mpc, pr);
   mpc_set_f(mpc, a, default_rnd);
}

inline mpcomplex::mpcomplex(const double a, mp_prec_t pr, mp_prec_t pi, mpc_rnd_t mode)
{
   mpc_init3(mpc, pr, pi);
   mpc_set_d(mpc, a, default_rnd);
}

inline mpcomplex::~mpcomplex() 
{ 
   mpc_clear(mpc);
}                           

inline mpcomplex& mpcomplex::operator=(const mpcomplex& a)
{
	if (this!= &a) mpc_set(mpc,a.mpc,default_rnd);
	return *this;
}

inline mpcomplex& mpcomplex::operator=(const mpc_t a)
{
	mpc_set(mpc,a,default_rnd);
	return *this;
}

inline mpcomplex& mpcomplex::operator=(const std::complex<double> a)
{
    mpc_set_d_d(mpc, a.real(), a.imag(), default_rnd);
    return *this;
}

inline mpcomplex& mpcomplex::operator=(const char* s)
{
     mpc_init3(mpc, default_real_prec, default_imag_prec);
     mpc_set_str(mpc, (char *)s,default_base,default_rnd);				
     return *this;
}

inline mpcomplex& mpcomplex::operator=(const mpreal& a)
{
   mpc_set_fr(mpc, (mpfr_ptr)(&a), default_rnd);
   return *this;
}

inline mpcomplex& mpcomplex::operator=(const mpfr_t a)
{
   mpc_set_fr(mpc, a, default_rnd);
   return *this;
}

inline mpcomplex& mpcomplex::operator=(const double a)
{
   mpc_set_d(mpc, a, default_rnd);
   return *this;
}
// + Addition
inline mpcomplex& mpcomplex::operator+=(const mpcomplex& a)
{
	mpc_add(mpc,mpc,a.mpc,default_rnd);
	return *this;
}

inline mpcomplex& mpcomplex::operator+=(const mpc_t a)
{
	*this += mpcomplex(a);
	return *this;
}

inline mpcomplex& mpcomplex::operator+=(const std::complex<double> a)
{
	return *this += mpcomplex(a);	
}

inline mpcomplex& mpcomplex::operator+=(const mpreal& a)
{
     mpc_add_fr(mpc, mpc, (mpfr_ptr)(&a), default_rnd);
     return *this;
}

inline mpcomplex& mpcomplex::operator+=(const mpfr_t a)
{
     mpc_add_fr(mpc, mpc, a, default_rnd);
     return *this;
}


inline mpcomplex& mpcomplex::operator+=(const double a)
{
     mpreal p = a;
     mpc_add_fr(mpc, mpc, (mpfr_ptr)(&p), default_rnd);
     return *this;
}

inline const mpcomplex mpcomplex::operator+()const
{
	return mpcomplex(*this);
}

inline mpcomplex& mpcomplex::operator++() 
{
	*this += 1.0;
	return *this;
}

inline const mpcomplex mpcomplex::operator++ (int)
{
	mpcomplex x(*this);
	*this += 1.0;
	return x;
}

// + Addition again
inline const mpcomplex operator+(const mpcomplex& a, const mpcomplex& b)
{
  if (!a.get_prec() == 0 && !b.get_prec() == 0) {
    if(a.get_prec()>b.get_prec()) return mpcomplex(a) += b;
    else return mpcomplex(b) += a;
  } else {
    mpcomplex tmp(a);
    tmp.set_prec2(std::max(a.get_prec_re(), b.get_prec_re()), std::max(a.get_prec_im(), b.get_prec_im()));
    return tmp += b;
  }
}

inline const mpcomplex operator+(const mpcomplex& a, const std::complex<double> b)
{
  return mpcomplex(a) += b;
}

inline const mpcomplex operator+(const mpcomplex& a, const char *b)
{
  return mpcomplex(b) += a;
}

inline const mpcomplex operator+(const std::complex<double> a, const mpcomplex& b)
{
  return mpcomplex(b) += a;
}

inline const mpcomplex operator+(const char *a, const mpcomplex& b)
{
  return mpcomplex(a) += b;
}


inline const mpcomplex operator+(const mpcomplex& a, const mpreal b)
{
  mpcomplex tmp(a);
  tmp.set_prec2(std::max(a.get_prec_re(), b.get_prec()), a.get_prec_im());
  return tmp += b;
}

inline const mpcomplex operator+(const mpreal a, const mpcomplex& b)
{
  mpcomplex tmp(b);
  tmp.set_prec2(std::max(b.get_prec_re(), a.get_prec()), b.get_prec_im());
  return tmp += a;
}

inline const mpcomplex operator+(const double a, const mpcomplex& b)
{
  mpcomplex tmp(b);
  return tmp += a;
}

inline const mpcomplex operator+(const int a, const mpcomplex& b)
{
  mpcomplex tmp(b);
  return tmp += a;
}

// - Subtraction
inline mpcomplex& mpcomplex::operator-=(const mpcomplex& a)
{
        mpc_sub(mpc,mpc,a.mpc,default_rnd);
        return *this;
}

inline mpcomplex& mpcomplex::operator-=(const mpc_t a)
{
	*this -= mpcomplex(a);
	return *this;
}

inline mpcomplex& mpcomplex::operator-=(const std::complex<double> a)
{
	return *this -= mpcomplex(a);	
}

inline mpcomplex& mpcomplex::operator-=(const mpreal& a)
{
     mpc_sub_fr(mpc, mpc, (mpfr_ptr)(&a), default_rnd);
     return *this;
}

inline mpcomplex& mpcomplex::operator-=(const mpfr_t a)
{
     mpc_sub_fr(mpc, mpc, a, default_rnd);
     return *this;
}


inline mpcomplex& mpcomplex::operator-=(const double a)
{
     mpreal p = a;
     mpc_sub_fr(mpc, mpc, (mpfr_ptr)(&p), default_rnd);
     return *this;
}

inline const mpcomplex mpcomplex::operator-()const
{
     mpcomplex u(*this);
     mpc_neg (u.mpc, u.mpc, default_rnd);
     return u;
}

inline mpcomplex& mpcomplex::operator--() 
{
	*this -= 1.0;
	return *this;
}

inline const mpcomplex mpcomplex::operator-- (int)
{
	mpcomplex x(*this);
	*this -= 1.0;
	return x;
}

inline const mpcomplex operator-(const mpcomplex& a, const mpcomplex& b)
{
  if (!a.get_prec() == 0 && !b.get_prec() == 0) {
    if(a.get_prec()>b.get_prec()) return mpcomplex(a) -= b;
    else return -(mpcomplex(b) -= a);
  } else {
    mpcomplex tmp(a);
    tmp.set_prec2(std::max(a.get_prec_re(), b.get_prec_re()), std::max(a.get_prec_im(), b.get_prec_im()));
    return tmp -= b;
  }
}

inline const mpcomplex operator-(const mpcomplex& a, const std::complex<double>& b)
{
    return -(mpcomplex(b) -= a);
}

inline const mpcomplex operator-(const mpcomplex& a, const char* b)
{
       return a-mpcomplex(b);
}

inline const mpcomplex operator-(const std::complex<double>& a, const mpcomplex& b)
{
    return mpcomplex(a) -= b;
}

inline const mpcomplex operator-(const char* a, const mpcomplex& b)
{
        return mpcomplex(a)-b;
}

inline const mpcomplex operator-(const mpcomplex& a, const mpreal b)
{
  if (!a.get_prec() == 0 && !b.get_prec() == 0) {
    if(a.get_prec()>b.get_prec()) return mpcomplex(a) -= b;
    else return -(mpcomplex(b) -= a);
  } else {
    mpcomplex tmp(a);
    tmp.set_prec2(std::max(a.get_prec_re(), b.get_prec()), a.get_prec_im());
    return tmp -= b;
  }
}

inline const mpcomplex operator-(const mpcomplex& a, const double b)
{
        return mpcomplex(a)-mpreal(b);
}

inline const mpcomplex operator-(const mpcomplex& a, const int b)
{
        return mpcomplex(a)-mpreal(b);
}

inline const mpcomplex operator-(const mpreal a, const mpcomplex& b)
{
        return mpcomplex(a)-mpcomplex(b);
}

inline const mpcomplex operator-(const double a, const mpcomplex& b)
{
        return mpcomplex(a)-mpcomplex(b);
}

inline const mpcomplex operator-(const int a, const mpcomplex& b)
{
        return mpcomplex(a)-mpcomplex(b);
}

// * multiplication
inline mpcomplex& mpcomplex::operator*=(const mpcomplex& a)
{
	mpc_mul(mpc,mpc,a.mpc,default_rnd);
	return *this;
}

inline mpcomplex& mpcomplex::operator*=(const mpc_t a)
{
	*this *= mpcomplex(a);
	return *this;
}

inline mpcomplex& mpcomplex::operator*=(const std::complex<double> a)
{
	return *this *= mpcomplex(a);	
}

inline mpcomplex& mpcomplex::operator*=(const mpreal& a)
{
     mpc_mul_fr(mpc, mpc, (mpfr_ptr)(&a), default_rnd);
     return *this;
}

inline mpcomplex& mpcomplex::operator*=(const mpfr_t a)
{
     mpc_mul_fr(mpc, mpc, a, default_rnd);
     return *this;
}


inline mpcomplex& mpcomplex::operator*=(const double a)
{
     mpreal p(a);
     mpc_add_fr(mpc, mpc, (mpfr_ptr)(&p), default_rnd);
     return *this;
}

inline const mpcomplex operator*(const mpcomplex& a, const mpcomplex& b)
{
  if (!a.get_prec() == 0 && !b.get_prec() == 0) {
    if(a.get_prec()>b.get_prec()) return mpcomplex(a) *= b;
    else return mpcomplex(b) *= a;
  } else {
    mpcomplex tmp(a);
    tmp.set_prec2(std::max(a.get_prec_re(), b.get_prec_re()), std::max(a.get_prec_im(), b.get_prec_im()));
    return tmp *= b;
  }
}

inline const mpcomplex operator*(const mpcomplex& a, const mpreal& b)
{
   mpcomplex tmp(a);
   tmp.set_prec2(std::max(a.get_prec_re(), b.get_prec()), a.get_prec_im());
   return tmp *= mpcomplex(b);
}

inline const mpcomplex operator*(const mpreal& a, const mpcomplex& b)
{
   mpcomplex tmp(a);
   tmp.set_prec2(std::max(b.get_prec_re(), a.get_prec()), b.get_prec_im());
   return tmp *= b;
}

// / division
inline mpcomplex& mpcomplex::operator/=(const mpcomplex& a)
{
	mpc_div(mpc,mpc,a.mpc,default_rnd);
	return *this;
}

inline mpcomplex& mpcomplex::operator/=(const mpc_t a)
{
	*this /= mpcomplex(a);
	return *this;
}

inline mpcomplex& mpcomplex::operator/=(const std::complex<double> a)
{
	return *this /= mpcomplex(a);	
}

inline mpcomplex& mpcomplex::operator/=(const mpreal& a)
{
     mpc_div_fr(mpc, mpc, (mpfr_ptr)(&a), default_rnd);
     return *this;
}

inline mpcomplex& mpcomplex::operator/=(const mpfr_t a)
{
     mpc_div_fr(mpc, mpc, a, default_rnd);
     return *this;
}

inline mpcomplex& mpcomplex::operator/=(const double a)
{
     mpreal p(a);
     mpc_div_fr(mpc, mpc, (mpfr_ptr)(&p), default_rnd);
     return *this;
}

inline bool operator == (const mpcomplex& a, const mpcomplex& b)
{
	return (mpc_cmp(a.mpc,b.mpc)==0);
}

inline bool operator == (const mpcomplex& a, const mpreal& b)
{
        mpcomplex c(b);
	return (mpc_cmp(a.mpc, c.mpc)==0);
}

inline bool operator == (const mpreal& a, const mpcomplex& b)
{
        mpcomplex c(a);
	return (mpc_cmp(c.mpc, b.mpc)==0);
}

inline bool operator != (const mpcomplex& a, const mpcomplex& b)
{
	return (!mpc_cmp(a.mpc, b.mpc)==0);
}

inline bool operator != (const mpcomplex& a, const mpreal& b)
{
        mpcomplex c(b);
	return (!mpc_cmp(a.mpc, c.mpc)==0);
}

inline bool operator != (const mpreal& a, const mpcomplex& b)
{
        mpcomplex c(a);
	return (!mpc_cmp(b.mpc, c.mpc)==0);
}

inline mpcomplex::operator std::complex<double>() const
{
  mpreal re, im;
  std::complex<double> tmp;
  re = (*this).real();
  im = (*this).imag();
  tmp.real() = mpfr_get_d(re,MPC_RND_RE(default_rnd));
  tmp.imag() = mpfr_get_d(im,MPC_RND_IM(default_rnd));
  return tmp;
}

inline const mpcomplex operator/(const mpcomplex& a, const mpcomplex& b)
{
  mpcomplex tmp(a);
  tmp.set_prec2(std::max(a.get_prec_re(), b.get_prec_re()), std::max(a.get_prec_im(), b.get_prec_im()));
  return tmp /= b;
}

inline const mpreal abs(const mpcomplex& a, mpfr_rnd_t rnd_mode)
{
        mpreal x;
        mpcomplex y(a);
	mpc_abs((mpfr_ptr)(&x), y.mpc, rnd_mode);
	return x;
}

inline const mpreal norm(const mpcomplex& a, mpfr_rnd_t rnd_mode)
{
        mpreal x;
        mpcomplex y(a);
	mpc_norm((mpfr_ptr)(&x), y.mpc, rnd_mode);
	return x;
}

inline const mpcomplex conj(const mpcomplex& a, mpc_rnd_t rnd_mode)
{
        mpcomplex y(a);
	mpc_conj(y.mpc, y.mpc, rnd_mode);
	return y;
}

inline const mpreal arg(const mpcomplex& a, mpfr_rnd_t rnd_mode)
{
        mpreal x(a.real());
        mpcomplex y(a);
	mpc_arg((mpfr_ptr)(&x), y.mpc, rnd_mode);
	return x;
}

inline const mpcomplex proj(const mpcomplex& a, mpc_rnd_t rnd_mode)
{
        mpcomplex y(a);
	mpc_proj(y.mpc, y.mpc, rnd_mode);
	return y;
}

inline const mpcomplex sqrt(const mpcomplex& a, mpc_rnd_t rnd_mode)
{
	mpcomplex x(a);
	mpc_sqrt(x.mpc, x.mpc, rnd_mode);
	return x;
}

inline const mpcomplex pow(const mpcomplex& a, const mpcomplex& b, mpc_rnd_t rnd_mode)
{
	mpcomplex x(a);
	mpc_pow(x.mpc, x.mpc, b.mpc, rnd_mode);
	return x;
}

inline const mpcomplex exp(const mpcomplex& a, mpc_rnd_t rnd_mode)
{
	mpcomplex x(a);
	mpc_exp(x.mpc, a.mpc, rnd_mode);
	return x;
}

inline const mpcomplex log(const mpcomplex& a, mpc_rnd_t rnd_mode)
{
	mpcomplex x(a);
	mpc_log(x.mpc, a.mpc, rnd_mode);
	return x;
}

inline const mpcomplex sin(const mpcomplex& a, mpc_rnd_t rnd_mode)
{
	mpcomplex x(a);
	mpc_sin(x.mpc, a.mpc, rnd_mode);
	return x;
}

inline const mpcomplex cos(const mpcomplex& a, mpc_rnd_t rnd_mode)
{
	mpcomplex x(a);
	mpc_cos(x.mpc, a.mpc, rnd_mode);
	return x;
}

inline const mpcomplex tan(const mpcomplex& a, mpc_rnd_t rnd_mode)
{
	mpcomplex x(a);
	mpc_tan(x.mpc, a.mpc, rnd_mode);
	return x;
}

inline const mpcomplex sinh(const mpcomplex& a, mpc_rnd_t rnd_mode)
{
	mpcomplex x(a);
	mpc_sinh(x.mpc, a.mpc, rnd_mode);
	return x;
}

inline const mpcomplex cosh(const mpcomplex& a, mpc_rnd_t rnd_mode)
{
	mpcomplex x(a);
	mpc_cosh(x.mpc, a.mpc, rnd_mode);
	return x;
}

inline const mpcomplex tanh(const mpcomplex& a, mpc_rnd_t rnd_mode)
{
	mpcomplex x(a);
	mpc_tanh(x.mpc, a.mpc, rnd_mode);
	return x;
}

inline const mpcomplex asin(const mpcomplex& a, mpc_rnd_t rnd_mode)
{
	mpcomplex x(a);
	mpc_asin(x.mpc, a.mpc, rnd_mode);
	return x;
}

inline const mpcomplex acos(const mpcomplex& a, mpc_rnd_t rnd_mode)
{
	mpcomplex x(a);
	mpc_acos(x.mpc, a.mpc, rnd_mode);
	return x;
}

inline const mpcomplex atan(const mpcomplex& a, mpc_rnd_t rnd_mode)
{
	mpcomplex x(a);
	mpc_atan(x.mpc, a.mpc, rnd_mode);
	return x;
}

inline const mpcomplex asinh(const mpcomplex& a, mpc_rnd_t rnd_mode)
{
	mpcomplex x(a);
	mpc_asinh(x.mpc, a.mpc, rnd_mode);
	return x;
}

inline const mpcomplex acosh(const mpcomplex& a, mpc_rnd_t rnd_mode)
{
	mpcomplex x(a);
	mpc_acosh(x.mpc, a.mpc, rnd_mode);
	return x;
}

inline const mpcomplex atanh(const mpcomplex& a, mpc_rnd_t rnd_mode)
{
	mpcomplex x(a);
	mpc_atanh(x.mpc, a.mpc, rnd_mode);
	return x;
}
inline const mpcomplex urandom_c (gmp_randstate_t& state)
{
	mpcomplex x;
	mpc_urandom(x.mpc,state);
	return x;
}

#if defined ___MPACK_BUILD_WITH_GMP___
inline mpcomplex::mpcomplex(const mpc_class& a)
{
   mp_prec_t pr, pi;
   pr = mpf_get_prec(a.real().get_mpf_t());  pi = mpf_get_prec(a.imag().get_mpf_t());
   mpc_init3(mpc, pr, pi);
   mpc_set_f_f(mpc, a.real().get_mpf_t(), a.imag().get_mpf_t(), default_rnd);
}

inline const mpcomplex operator-(const mpcomplex& a, const mpc_class& b)
{
   mpcomplex tmp(b);
   return mpcomplex(a) -= b;
}

inline const mpcomplex operator-(const mpc_class& a, const mpcomplex& b)
{
   mpcomplex tmp(a);
   mpcomplex tmp2(b);
   return -(mpcomplex(a) -= b);
}

inline mpc_class cast2mpc_class(const mpcomplex &b)
{
//mpcomplex -> mpfr, mpfr -> mpf, mpf -> mpc_class
  mpf_t  tmpre, tmpim;
  mpfr_t tmpre2, tmpim2;
  mpcomplex a(b);
  mp_prec_t pr, pi;

  pr = a.get_prec_re();
  pi = a.get_prec_im();
  
  mpfr_init2(tmpre2, pr); mpfr_init2(tmpim2, pi);
  mpfr_set(tmpre2, (mpfr_ptr)a.real(), MPC_RND_RE(mpcomplex::default_rnd));
  mpfr_set(tmpim2, (mpfr_ptr)a.imag(), MPC_RND_IM(mpcomplex::default_rnd));

  mpf_init2 (tmpre, pr); mpf_init2 (tmpim, pr);
  mpfr_get_f(tmpre, tmpre2, MPC_RND_RE(mpcomplex::default_rnd));  mpfr_get_f(tmpim, tmpim2, MPC_RND_IM(mpcomplex::default_rnd));

  mpc_class tmp(tmpre, tmpim);

  mpf_clear(tmpre);
  mpf_clear(tmpim);

  mpfr_clear(tmpre2);
  mpfr_clear(tmpim2);
  return tmp;
}     

#endif

#if defined ___MPACK_BUILD_WITH_QD___
inline mpcomplex::mpcomplex(const qd_complex& a)
{
   mpfr_t mp_real, mp_imag;
   mpc_init3(mpc, default_real_prec, default_imag_prec);

   mpfr_init2(mp_real, default_real_prec);
   mpfr_set_d (mp_real, a.real().x[0],  MPC_RND_RE(default_rnd));
   mpfr_add_d (mp_real, mp_real, a.real().x[1], MPC_RND_RE(default_rnd));
   mpfr_add_d (mp_real, mp_real, a.real().x[2], MPC_RND_RE(default_rnd));
   mpfr_add_d (mp_real, mp_real, a.real().x[3], MPC_RND_RE(default_rnd));

   mpfr_init2(mp_imag, default_imag_prec);
   mpfr_set_d (mp_imag, a.imag().x[0], MPC_RND_IM(default_rnd));
   mpfr_add_d (mp_imag, mp_imag, a.imag().x[1], MPC_RND_IM(default_rnd));
   mpfr_add_d (mp_imag, mp_imag, a.imag().x[2], MPC_RND_IM(default_rnd));
   mpfr_add_d (mp_imag, mp_imag, a.imag().x[3], MPC_RND_IM(default_rnd));
   
   mpc_set_fr_fr(mpc, mp_real, mp_imag, default_rnd);
   mpfr_clear(mp_imag);
   mpfr_clear(mp_real);
}
inline qd_complex cast2qd_complex(const mpcomplex &b)
{
  std::complex<double> p[4];
  mpcomplex c (b);
  qd_complex q;

  p[0] = c;
  c = c - p[0];
  p[1] = c;
  c = c - p[1];
  p[2] = c;
  c = c - p[2];
  p[3] = c;
//do not optimize
  q = p[0];
  q = q + p[1];
  q = q + p[2];
  q = q + p[3];
  return q;
}     
inline const mpcomplex operator-(const mpcomplex& a, const qd_complex& b)
{
   return mpcomplex(b) -= a;
}

inline const mpcomplex operator-(const qd_complex& a, const mpcomplex& b)
{
   return -(mpcomplex(a) -= b);
}

inline mpcomplex& mpcomplex::operator=(const qd_complex& a)
{
   mpcomplex tmp(a);
   *this = tmp; 
   return *this;
}

#endif

#if defined ___MPACK_BUILD_WITH_DD___
inline mpcomplex::mpcomplex(const dd_complex& a)
{
   mpfr_t mp_real, mp_imag;
   mpc_init3(mpc, default_real_prec, default_imag_prec);

   mpfr_init2(mp_real, default_real_prec);
   mpfr_set_d (mp_real, a.real().x[0],  MPC_RND_RE(default_rnd));
   mpfr_add_d (mp_real, mp_real, a.real().x[1], MPC_RND_RE(default_rnd));

   mpfr_init2(mp_imag, default_imag_prec);
   mpfr_set_d (mp_imag, a.imag().x[0], MPC_RND_IM(default_rnd));
   mpfr_add_d (mp_imag, mp_imag, a.imag().x[1], MPC_RND_IM(default_rnd));
   
   mpc_set_fr_fr(mpc, mp_real, mp_imag, default_rnd);
   mpfr_clear(mp_imag);
   mpfr_clear(mp_real);
}
inline dd_complex cast2dd_complex(const mpcomplex &b)
{
  std::complex<double> p[2];
  mpcomplex c (b);
  dd_complex q;

  p[0] = c;
  c = c - p[0];
  p[1] = c;
  c = c - p[1];
//do not optimize
  q = p[0];
  q = q + p[1];
  return q;
}     
inline const mpcomplex operator-(const mpcomplex& a, const dd_complex& b)
{
   return mpcomplex(b) -= a;
}

inline const mpcomplex operator-(const dd_complex& a, const mpcomplex& b)
{
   return -(mpcomplex(a) -= b);
}

inline mpcomplex& mpcomplex::operator=(const dd_complex& a)
{
   mpcomplex tmp(a);
   *this = tmp; 
   return *this;
}

#endif

}

#endif /* __MP_COMPLEX_H__ */

[Mplapack-devel] Bug in Rgerfs

[Fletcher, J. P <j.p...@as...>]

// solve_dd_refine.cpp from  solve_double_refine.cpp
// solve_dd.cpp from inv_mat_dd.cpp
//
// Get inverse of Matrix A via Rgetri, using double
// This file is freely usable.
// written by Nakata Maho, 2009/9/23.

#include <mblas_dd.h>
#include <mlapack_dd.h>
#include <stdio.h>

//Matlab/Octave format
void printmat(int N, int M, dd_real *A, int LDA)
{
    dd_real mtmp;

    printf("[ ");
    for (int i = 0; i < N; i++) {
	printf("[ ");
	for (int j = 0; j < M; j++) {
	    mtmp = A[i + j * LDA];
	    std::cout << mtmp;
	    //	    printf("%5.2e", mtmp);
	    if (j < M - 1)
		printf(", ");
	}
	if (i < N - 1)
	    printf("]; ");
	else
	    printf("] ");
    }
    printf("]");
}

int main()
{
    mpackint n = 3;
    mpackint lwork, info;

    dd_real *A = new dd_real[n * n];
    dd_real *AF = new dd_real[n * n];
    mpackint *ipiv = new mpackint[n];
    mpackint *ipiv1 = new mpackint[n+1];
    mpackint *iwork = new mpackint[n+1];

    mpackint nrhs = 1;
 
    dd_real *B = new dd_real[n*nrhs];
    dd_real *X = new dd_real[n*nrhs];
    dd_real *B1 = new dd_real[n*nrhs+1];
    dd_real *X1 = new dd_real[n*nrhs+1];
    dd_real *work = new dd_real[3 * n + 1];
    dd_real *berr = new dd_real[n];
    dd_real *ferr = new dd_real[n];
//setting A matrix
    A[0 + 0 * n] = 1;    A[0 + 1 * n] = 4;    A[0 + 2 * n] = 6;
    A[1 + 0 * n] = 2;    A[1 + 1 * n] = 9;    A[1 + 2 * n] = 14;
    A[2 + 0 * n] = 3;    A[2 + 1 * n] = 14;   A[2 + 2 * n] = 23;
 
    B[0] = 1; B[1] = 2; B[2] = 3;

    printf("A = ");
    printmat(n, n, A, n);
    printf("\n");
    printf("B = ");
    printmat(n, nrhs, B, n);
    printf("\n");

    for (int i = 0; i < n*n; i++) AF[i] = A[i];
    for (int i = 0; i < n*nrhs; i++) X[i] = B[i];

    Rgetrf(n, n, AF, n, ipiv, &info);
    Rgetrs("n", n, nrhs, AF, n, ipiv, X, n, &info);

    printf("X = ");
    printmat(n, nrhs, X, n);
    printf("\n");
 
    for (int i = 0; i < n; i++) ipiv1[i+1] = ipiv[i];
    for (int i = 0; i < n; i++) B1[i+1] = B[i];
    for (int i = 0; i < n; i++) X1[i+1] = X[i];

    Rgerfs("n", n, nrhs, A, n, AF, n, ipiv1, B1, n, X1, n,
           ferr, berr, work, iwork, &info);
    printf("info = %d\n",info);
    printf("revised X = ");
    printmat(n, nrhs, X, n);
    printf("\n");
    printmat(nrhs, 1, ferr, n);
    printmat(nrhs, 1, berr, n);
    printf("\n");

    delete[]work;
    delete[]iwork;
    delete[]ipiv;
    delete[]ipiv1;
    delete[]A;
    delete[]AF;
    delete[]B;
    delete[]X;
    delete[]berr;
    delete[]ferr;


}

Re: [Mplapack-devel] FW: Fortran examples

[Maho N. <ch...@ma...>]

From: Maho NAKATA <ch...@ma...>
Subject: Re: [Mplapack-devel] FW: Fortran examples
Date: Tue, 05 Apr 2011 11:08:52 +0900 (JST)

> Hi,
> 
> From: "Fletcher, John P" <j.p...@as...>
> Subject: RE: [Mplapack-devel] FW: Fortran examples
> Date: Fri, 01 Apr 2011 09:58:54 +0100
> 
>> Nakata Maho
>> 
>> I was looking at the choice last night.  I think I am going to use mpfr (mpreal C++ interface) and mpcomplex.  I needed to reorder the libraries, -lmpc before -lmpfr as a test case was not linking properly.  I had copied MPFRLIBS in the examples/mlapack Makefile.
> 
> thanks for your comment. Could you please tell me which part of MPACK is wrong?
> 
>> There does not seem to be an operator<< in mcomplex but it should not be much work for me to make one.
> 
> Could you please tell me what "operator <<" do for mpcomplex?
> 
>> I am using gfortran, so plan to follow the MAGMA style to make an interface, but I will start with some C++ examples with MPACK to get a better understanding. I am more familiar with LAPACK than MPFR.
> 
> ok. Please let me know when idea comes up. Gfortran interface is also very important.
> 
>> You will find I am quite active at the moment on the MAGMA forum.
>> 
>> Are you in touch with Goto Kazusige?  There is unfortunately a bug of some sort in the Nehalem version of GotoBlas2 which means that safe operation requires use of the 4 core version even on an 8 core machine. It does not look as though maintenance is being done any more on the distribution site.
>> 
>> I don't have full details but can get them if needed.
> 
> I e-mailed to Goto Kazushige-san to ask your question.


Goto san kindly told me that even though logical core is counted as eight, 
GotoBLAS2 enumerate the number of physical core and set it as the number of threads.

thanks
-- Nakata Maho http://accc.riken.jp/maho/ , JA OOO http://ja.openoffice.org/
http://blog.goo.ne.jp/nakatamaho/ ,GPG: http://accc.riken.jp/maho/maho.pgp.txt

Re: [Mplapack-devel] FW: Fortran examples

[Maho N. <ch...@ma...>]

Hi,

From: "Fletcher, John P" <j.p...@as...>
Subject: RE: [Mplapack-devel] FW: Fortran examples
Date: Fri, 01 Apr 2011 09:58:54 +0100

> Nakata Maho
> 
> I was looking at the choice last night.  I think I am going to use mpfr (mpreal C++ interface) and mpcomplex.  I needed to reorder the libraries, -lmpc before -lmpfr as a test case was not linking properly.  I had copied MPFRLIBS in the examples/mlapack Makefile.

thanks for your comment. Could you please tell me which part of MPACK is wrong?

> There does not seem to be an operator<< in mcomplex but it should not be much work for me to make one.

Could you please tell me what "operator <<" do for mpcomplex?

> I am using gfortran, so plan to follow the MAGMA style to make an interface, but I will start with some C++ examples with MPACK to get a better understanding. I am more familiar with LAPACK than MPFR.

ok. Please let me know when idea comes up. Gfortran interface is also very important.

> You will find I am quite active at the moment on the MAGMA forum.
> 
> Are you in touch with Goto Kazusige?  There is unfortunately a bug of some sort in the Nehalem version of GotoBlas2 which means that safe operation requires use of the 4 core version even on an 8 core machine. It does not look as though maintenance is being done any more on the distribution site.
> 
> I don't have full details but can get them if needed.

I e-mailed to Goto Kazushige-san to ask your question.

Thanks


> Dr John P. Fletcher Tel: (44) 121 204 3389 (direct line), FAX: (44) 121 204 3678
> Chemical Engineering and Applied Chemistry (CEAC),
> Associate Dean - External Relations,
> School of Engineering and Applied Science (EAS),
> Aston University, Aston Triangle, BIRMINGHAM B4 7ET  U.K.               
> 
> 
> 
> -----Original Message-----
> From: Maho NAKATA [mailto:mah...@gm...] On Behalf Of Maho NAKATA
> Sent: 01 April 2011 01:32
> To: Fletcher, John P
> Cc: mpl...@li...
> Subject: Re: [Mplapack-devel] FW: Fortran examples
> 
> Hi Fletcher,
> 
> Oh you're developing MAGMA! Great. It is very fast on C2050.
> 
> Anyway, interface to FORTRAN77 would be impossible, but FORTRAN90 or 95 might be
> possible. Which MP library do you want to use? GMP, MPFR, QD?
> thanks,
>  Nakata Maho
> 
> From: "Fletcher, John P" <j.p...@as...>
> Subject: [Mplapack-devel] FW: Fortran examples
> Date: Thu, 31 Mar 2011 13:05:04 +0100
> 
>> Maho
>> 
>> On thinking about it, it was the wrong question to ask as the structure is such that C++ types will not have equivalents in FORTRAN.  The problem is that I want to adapt some code which is in FORTRAN so I am going to have to find some way to interface the two.  I am working with MAGMA (LAPACK for GPGPUs see http://icl.cs.utk.edu/magma/index.html) and I have to interface FORTRAN to C++ there.
>> 
>> Thank you
>> 
>> John
>> 
>> Dr John P. Fletcher Tel: (44) 121 204 3389 (direct line), FAX: (44) 121 204 3678
>> Chemical Engineering and Applied Chemistry (CEAC),
>> Associate Dean - External Relations,
>> School of Engineering and Applied Science (EAS),
>> Aston University, Aston Triangle, BIRMINGHAM B4 7ET  U.K.               
>> 
>> 
>> 
>> -----Original Message-----
>> From: Maho NAKATA [mailto:ch...@ma...] 
>> Sent: 31 March 2011 03:33
>> To: Fletcher, John P
>> Cc: mpl...@li...
>> Subject: Re: [Mplapack-devel] Fortran examples
>> 
>> Hi Fletcher, 
>> 
>> From: "Fletcher, John P" <j.p...@as...>
>> Subject: [Mplapack-devel] Fortran examples
>> Date: Wed, 30 Mar 2011 14:07:36 +0100
>> 
>>> I have just started using MPACK and successfully installed it on Ubuntu Linux 10.4 (64 bit).
>> 
>> Congratulations!
>> 
>>> Please, are there any examples in FORTRAN available?
>> 
>> Unfortunately, no...
>> 
>> thanks
>> -- Nakata Maho http://accc.riken.jp/maho/ , JA OOO http://ja.openoffice.org/
>> http://blog.goo.ne.jp/nakatamaho/ ,GPG: http://accc.riken.jp/maho/maho.pgp.txt
>> 
>> ------------------------------------------------------------------------------
>> Create and publish websites with WebMatrix
>> Use the most popular FREE web apps or write code yourself; 
>> WebMatrix provides all the features you need to develop and 
>> publish your website. http://p.sf.net/sfu/ms-webmatrix-sf
>> _______________________________________________
>> Mplapack-devel mailing list
>> Mpl...@li...
>> https://lists.sourceforge.net/lists/listinfo/mplapack-devel
>> 
>> ------------------------------------------------------------------------------
>> Create and publish websites with WebMatrix
>> Use the most popular FREE web apps or write code yourself; 
>> WebMatrix provides all the features you need to develop and 
>> publish your website. http://p.sf.net/sfu/ms-webmatrix-sf
>> _______________________________________________
>> Mplapack-devel mailing list
>> Mpl...@li...
>> https://lists.sourceforge.net/lists/listinfo/mplapack-devel
>> 
> 

[Mplapack-devel] Versions of MPFR and GMP on Ubuntu Linux

[Fletcher, J. P <j.p...@as...>]

I have tried out dynamic linking and find that I get warnings associated with the versions of mpfr and gmp libraries which are being used.

g++ -c -O2 -fopenmp -I/home/fletcher/MPACK/include -I/home/fletcher/MPACK/include/mpack -I/home/fletcher/MPACK/include/qd trial_mpcomplex.cpp
g++ -o  trial_mpcomplex trial_mpcomplex.o -fopenmp -L/home/fletcher/MPACK/lib -lmlapack_mpfr -lmblas_mpfr -lmpfrcxx -lmpc -lmpfr -lgmp
/usr/bin/ld: warning: libmpfr.so.4, needed by /usr/local/lib/libmpc.so, may conflict with libmpfr.so.1
/usr/bin/ld: warning: libgmp.so.10, needed by /usr/local/lib/libmpc.so, may conflict with libgmp.so.3
On inspection libmpr.so.1 and libgmp.so.3 are located in /usr/lib and are the ones installed by the operating system:

fletcher@fletcher-desktop:~$<mailto:fletcher@fletcher-desktop:~$> apt-cache policy libgmp3c2
libgmp3c2:
  Installed: 2:4.3.2+dfsg-1ubuntu1
  Candidate: 2:4.3.2+dfsg-1ubuntu1
  Version table:
 *** 2:4.3.2+dfsg-1ubuntu1 0
        500 http://gb.archive.ubuntu.com/ubuntu/ lucid/main Packages
        100 /var/lib/dpkg/status
fletcher@fletcher-desktop:~$<mailto:fletcher@fletcher-desktop:~$> apt-cache policy libmpfr1ldbl
libmpfr1ldbl:
  Installed: 2.4.2-3ubuntu1
  Candidate: 2.4.2-3ubuntu1
  Version table:
 *** 2.4.2-3ubuntu1 0
        500 http://gb.archive.ubuntu.com/ubuntu/ lucid/main Packages
        100 /var/lib/dpkg/status
The ones in /usr/local/lib I think have been installed by MPACK (libgmp.so.10 and libmpfr.so.4) as they were updated when I did a rebuild of MPACK.

It is not unusual to have this sort of problem with Ubuntu as it follows Debian in using /usr/lib where other systems use /usr/local/lib.

Is there anything I can do to resolve this.

I don't want to remove the versions in /usr/lib as they are likely to be used by other software, such as gcc.

Thank you

John

Re: [Mplapack-devel] FW: Fortran examples

[Fletcher, J. P <j.p...@as...>]

Nakata Maho

I was looking at the choice last night.  I think I am going to use mpfr (mpreal C++ interface) and mpcomplex.  I needed to reorder the libraries, -lmpc before -lmpfr as a test case was not linking properly.  I had copied MPFRLIBS in the examples/mlapack Makefile.

There does not seem to be an operator<< in mcomplex but it should not be much work for me to make one.

I am using gfortran, so plan to follow the MAGMA style to make an interface, but I will start with some C++ examples with MPACK to get a better understanding. I am more familiar with LAPACK than MPFR.

You will find I am quite active at the moment on the MAGMA forum.

Are you in touch with Goto Kazusige?  There is unfortunately a bug of some sort in the Nehalem version of GotoBlas2 which means that safe operation requires use of the 4 core version even on an 8 core machine. It does not look as though maintenance is being done any more on the distribution site.

I don't have full details but can get them if needed.

Best wishes

John Fletcher

Dr John P. Fletcher Tel: (44) 121 204 3389 (direct line), FAX: (44) 121 204 3678
Chemical Engineering and Applied Chemistry (CEAC),
Associate Dean - External Relations,
School of Engineering and Applied Science (EAS),
Aston University, Aston Triangle, BIRMINGHAM B4 7ET  U.K.               



-----Original Message-----
From: Maho NAKATA [mailto:mah...@gm...] On Behalf Of Maho NAKATA
Sent: 01 April 2011 01:32
To: Fletcher, John P
Cc: mpl...@li...
Subject: Re: [Mplapack-devel] FW: Fortran examples

Hi Fletcher,

Oh you're developing MAGMA! Great. It is very fast on C2050.

Anyway, interface to FORTRAN77 would be impossible, but FORTRAN90 or 95 might be
possible. Which MP library do you want to use? GMP, MPFR, QD?
thanks,
 Nakata Maho

From: "Fletcher, John P" <j.p...@as...>
Subject: [Mplapack-devel] FW: Fortran examples
Date: Thu, 31 Mar 2011 13:05:04 +0100

> Maho
> 
> On thinking about it, it was the wrong question to ask as the structure is such that C++ types will not have equivalents in FORTRAN.  The problem is that I want to adapt some code which is in FORTRAN so I am going to have to find some way to interface the two.  I am working with MAGMA (LAPACK for GPGPUs see http://icl.cs.utk.edu/magma/index.html) and I have to interface FORTRAN to C++ there.
> 
> Thank you
> 
> John
> 
> Dr John P. Fletcher Tel: (44) 121 204 3389 (direct line), FAX: (44) 121 204 3678
> Chemical Engineering and Applied Chemistry (CEAC),
> Associate Dean - External Relations,
> School of Engineering and Applied Science (EAS),
> Aston University, Aston Triangle, BIRMINGHAM B4 7ET  U.K.               
> 
> 
> 
> -----Original Message-----
> From: Maho NAKATA [mailto:ch...@ma...] 
> Sent: 31 March 2011 03:33
> To: Fletcher, John P
> Cc: mpl...@li...
> Subject: Re: [Mplapack-devel] Fortran examples
> 
> Hi Fletcher, 
> 
> From: "Fletcher, John P" <j.p...@as...>
> Subject: [Mplapack-devel] Fortran examples
> Date: Wed, 30 Mar 2011 14:07:36 +0100
> 
>> I have just started using MPACK and successfully installed it on Ubuntu Linux 10.4 (64 bit).
> 
> Congratulations!
> 
>> Please, are there any examples in FORTRAN available?
> 
> Unfortunately, no...
> 
> thanks
> -- Nakata Maho http://accc.riken.jp/maho/ , JA OOO http://ja.openoffice.org/
> http://blog.goo.ne.jp/nakatamaho/ ,GPG: http://accc.riken.jp/maho/maho.pgp.txt
> 
> ------------------------------------------------------------------------------
> Create and publish websites with WebMatrix
> Use the most popular FREE web apps or write code yourself; 
> WebMatrix provides all the features you need to develop and 
> publish your website. http://p.sf.net/sfu/ms-webmatrix-sf
> _______________________________________________
> Mplapack-devel mailing list
> Mpl...@li...
> https://lists.sourceforge.net/lists/listinfo/mplapack-devel
> 
> ------------------------------------------------------------------------------
> Create and publish websites with WebMatrix
> Use the most popular FREE web apps or write code yourself; 
> WebMatrix provides all the features you need to develop and 
> publish your website. http://p.sf.net/sfu/ms-webmatrix-sf
> _______________________________________________
> Mplapack-devel mailing list
> Mpl...@li...
> https://lists.sourceforge.net/lists/listinfo/mplapack-devel
> 

Re: [Mplapack-devel] FW: Fortran examples

[Maho N. <ch...@ma...>]

Hi Fletcher,

Oh you're developing MAGMA! Great. It is very fast on C2050.

Anyway, interface to FORTRAN77 would be impossible, but FORTRAN90 or 95 might be
possible. Which MP library do you want to use? GMP, MPFR, QD?
thanks,
 Nakata Maho

From: "Fletcher, John P" <j.p...@as...>
Subject: [Mplapack-devel] FW: Fortran examples
Date: Thu, 31 Mar 2011 13:05:04 +0100

> Maho
> 
> On thinking about it, it was the wrong question to ask as the structure is such that C++ types will not have equivalents in FORTRAN.  The problem is that I want to adapt some code which is in FORTRAN so I am going to have to find some way to interface the two.  I am working with MAGMA (LAPACK for GPGPUs see http://icl.cs.utk.edu/magma/index.html) and I have to interface FORTRAN to C++ there.
> 
> Thank you
> 
> John
> 
> Dr John P. Fletcher Tel: (44) 121 204 3389 (direct line), FAX: (44) 121 204 3678
> Chemical Engineering and Applied Chemistry (CEAC),
> Associate Dean - External Relations,
> School of Engineering and Applied Science (EAS),
> Aston University, Aston Triangle, BIRMINGHAM B4 7ET  U.K.               
> 
> 
> 
> -----Original Message-----
> From: Maho NAKATA [mailto:ch...@ma...] 
> Sent: 31 March 2011 03:33
> To: Fletcher, John P
> Cc: mpl...@li...
> Subject: Re: [Mplapack-devel] Fortran examples
> 
> Hi Fletcher, 
> 
> From: "Fletcher, John P" <j.p...@as...>
> Subject: [Mplapack-devel] Fortran examples
> Date: Wed, 30 Mar 2011 14:07:36 +0100
> 
>> I have just started using MPACK and successfully installed it on Ubuntu Linux 10.4 (64 bit).
> 
> Congratulations!
> 
>> Please, are there any examples in FORTRAN available?
> 
> Unfortunately, no...
> 
> thanks
> -- Nakata Maho http://accc.riken.jp/maho/ , JA OOO http://ja.openoffice.org/
> http://blog.goo.ne.jp/nakatamaho/ ,GPG: http://accc.riken.jp/maho/maho.pgp.txt
> 
> ------------------------------------------------------------------------------
> Create and publish websites with WebMatrix
> Use the most popular FREE web apps or write code yourself; 
> WebMatrix provides all the features you need to develop and 
> publish your website. http://p.sf.net/sfu/ms-webmatrix-sf
> _______________________________________________
> Mplapack-devel mailing list
> Mpl...@li...
> https://lists.sourceforge.net/lists/listinfo/mplapack-devel
> 
> ------------------------------------------------------------------------------
> Create and publish websites with WebMatrix
> Use the most popular FREE web apps or write code yourself; 
> WebMatrix provides all the features you need to develop and 
> publish your website. http://p.sf.net/sfu/ms-webmatrix-sf
> _______________________________________________
> Mplapack-devel mailing list
> Mpl...@li...
> https://lists.sourceforge.net/lists/listinfo/mplapack-devel
> 

[Mplapack-devel] FW: Fortran examples

[Fletcher, J. P <j.p...@as...>]

Maho

On thinking about it, it was the wrong question to ask as the structure is such that C++ types will not have equivalents in FORTRAN.  The problem is that I want to adapt some code which is in FORTRAN so I am going to have to find some way to interface the two.  I am working with MAGMA (LAPACK for GPGPUs see http://icl.cs.utk.edu/magma/index.html) and I have to interface FORTRAN to C++ there.

Thank you

John

Dr John P. Fletcher Tel: (44) 121 204 3389 (direct line), FAX: (44) 121 204 3678
Chemical Engineering and Applied Chemistry (CEAC),
Associate Dean - External Relations,
School of Engineering and Applied Science (EAS),
Aston University, Aston Triangle, BIRMINGHAM B4 7ET  U.K.               



-----Original Message-----
From: Maho NAKATA [mailto:ch...@ma...] 
Sent: 31 March 2011 03:33
To: Fletcher, John P
Cc: mpl...@li...
Subject: Re: [Mplapack-devel] Fortran examples

Hi Fletcher, 

From: "Fletcher, John P" <j.p...@as...>
Subject: [Mplapack-devel] Fortran examples
Date: Wed, 30 Mar 2011 14:07:36 +0100

> I have just started using MPACK and successfully installed it on Ubuntu Linux 10.4 (64 bit).

Congratulations!

> Please, are there any examples in FORTRAN available?

Unfortunately, no...

thanks
-- Nakata Maho http://accc.riken.jp/maho/ , JA OOO http://ja.openoffice.org/
http://blog.goo.ne.jp/nakatamaho/ ,GPG: http://accc.riken.jp/maho/maho.pgp.txt

------------------------------------------------------------------------------
Create and publish websites with WebMatrix
Use the most popular FREE web apps or write code yourself; 
WebMatrix provides all the features you need to develop and 
publish your website. http://p.sf.net/sfu/ms-webmatrix-sf
_______________________________________________
Mplapack-devel mailing list
Mpl...@li...
https://lists.sourceforge.net/lists/listinfo/mplapack-devel

Re: [Mplapack-devel] Fortran examples

[Maho N. <ch...@ma...>]

Hi Fletcher, 

From: "Fletcher, John P" <j.p...@as...>
Subject: [Mplapack-devel] Fortran examples
Date: Wed, 30 Mar 2011 14:07:36 +0100

> I have just started using MPACK and successfully installed it on Ubuntu Linux 10.4 (64 bit).

Congratulations!

> Please, are there any examples in FORTRAN available?

Unfortunately, no...

thanks
-- Nakata Maho http://accc.riken.jp/maho/ , JA OOO http://ja.openoffice.org/
http://blog.goo.ne.jp/nakatamaho/ ,GPG: http://accc.riken.jp/maho/maho.pgp.txt

[Mplapack-devel] Fortran examples

[Fletcher, J. P <j.p...@as...>]

I have just started using MPACK and successfully installed it on Ubuntu Linux 10.4 (64 bit).

Please, are there any examples in FORTRAN available?

Thanks

John

Dr John P. Fletcher Tel: (44) 121 204 3389 (direct line), FAX: (44) 121 204 3678
Chemical Engineering and Applied Chemistry (CEAC),
Associate Dean - External Relations,
School of Engineering and Applied Science (EAS),
Aston University, Aston Triangle, BIRMINGHAM B4 7ET  U.K.