Leasqr matrix singular to machine precision

[Archambault Fabien]

Hi,

I am working on octave since a few days because I need to fit a 
potential on Lennard-Jones parameters (for those who understand).

So I created octave files that can be run in bash.

I think my equations are well written but if possible someone can check.

The script is :

#!/usr/bin/octave -qf
source "octave.script"
x = [ 1.367; 0.224500; 1.768200; -0.120; -0.046000; -0.152100 ];
r(1,1) = RTEMP(1,14);
r(1,2) = RTEMP(2,14);
r(1,3) = RTEMP(3,14);
y(1) = DELTA(14);
r(2,1) = RTEMP(1,1);
r(2,2) = RTEMP(2,1);
r(2,3) = RTEMP(3,1);
y(2) = DELTA(1);
r(3,1) = RTEMP(1,5);
r(3,2) = RTEMP(2,5);
r(3,3) = RTEMP(3,5);
y(3) = DELTA(5);
r(4,1) = RTEMP(1,13);
r(4,2) = RTEMP(2,13);
r(4,3) = RTEMP(3,13);
y(4) = DELTA(13);
r(5,1) = RTEMP(1,9);
r(5,2) = RTEMP(2,9);
r(5,3) = RTEMP(3,9);
y(5) = DELTA(9);
r(6,1) = RTEMP(1,12);
r(6,2) = RTEMP(2,12);
r(6,3) = RTEMP(3,12);
y(6) = DELTA(12);
function y = f(r,x)
    y = sqrt(x(4)*x(5)) * ( ((x(1)/2+x(2)/2) ./ r(:,1)).^12 - 2 * 
((x(1)/2+x(2)/2) ./ r(:,1)).^6 ) \
          + sqrt(x(4)*x(5)) * ( ((x(1)/2+x(2)/2) ./ r(:,2)).^12 - 2 * 
((x(1)/2+x(2)/2) ./ r(:,2)).^6 ) \
          + sqrt(x(4)*x(6)) * ( ((x(1)/2+x(3)/2) ./ r(:,3)).^12 - 2 * 
((x(1)/2+x(3)/2) ./ r(:,3)).^6 );
endfunction
[f1,result,kvg,iter] = leasqr(r,y,x,"f")

See in attach file for octave.script.

When I launch the script the message is :

warning: in /usr/share/octave/2.1.73/site/m/octave-forge/optim/leasqr.m 
near line 299, column 5:
warning: division by zero
warning: division by zero
warning: inverse: matrix singular to machine precision, rcond = 6.93366e-19
f1 =

  -4.8789e-04
  -7.2189e-02
  -2.5896e-02
  -9.8740e-04
  -1.0229e-02
  -2.1931e-03

result =

  1.367000
  0.224500
  1.768200
 -0.120000
 -0.046000
 -0.152100

kvg = 1
iter = 1

The version of octave is 2.1.73 on Mandriva 2007.1. The script leasqr is 
taken from octave-forge (in rpm files and tested with script available 
online).

I have searched also on the mailing list but did not find any clue for 
this problem.

Thanks in advance for your answer,
Fabien Archambault

--
Fabien Archambault
Equipe de dynamique des assemblages membranaires
Unité Mixte de Recherches CNRS UHP 7565
Université Henri-Poincaré, Nancy I BP 239,
54506 Vandoeuvre-lès-Nancy, cedex France

Tél : 03.83.68.43.96

GROUPANOM1 = [ 20 ];
GROUPAX1 = [ 0 ];
GROUPAY1 = [ 0 ];
GROUPAZ1 = [ 0 ];
GROUPANOM2 = [ 20 ];
GROUPAX2 = [ 0 ];
GROUPAY2 = [ 0 ];
GROUPAZ2 = [ 0 ];
GROUPANOM3 = [ 20 ];
GROUPAX3 = [ 0 ];
GROUPAY3 = [ 0 ];
GROUPAZ3 = [ 0 ];
GROUPANOM4 = [ 20 ];
GROUPAX4 = [ 0 ];
GROUPAY4 = [ 0 ];
GROUPAZ4 = [ 0 ];
GROUPANOM5 = [ 20 ];
GROUPAX5 = [ 0 ];
GROUPAY5 = [ 0 ];
GROUPAZ5 = [ 0 ];
GROUPANOM6 = [ 20 ];
GROUPAX6 = [ 0 ];
GROUPAY6 = [ 0 ];
GROUPAZ6 = [ 0 ];
GROUPANOM7 = [ 20 ];
GROUPAX7 = [ 0 ];
GROUPAY7 = [ 0 ];
GROUPAZ7 = [ 0 ];
GROUPANOM8 = [ 20 ];
GROUPAX8 = [ 0 ];
GROUPAY8 = [ 0 ];
GROUPAZ8 = [ 0 ];
GROUPANOM9 = [ 20 ];
GROUPAX9 = [ 0 ];
GROUPAY9 = [ 0 ];
GROUPAZ9 = [ 0 ];
GROUPANOM10 = [ 20 ];
GROUPAX10 = [ 0 ];
GROUPAY10 = [ 0 ];
GROUPAZ10 = [ 0 ];
GROUPANOM11 = [ 20 ];
GROUPAX11 = [ 0 ];
GROUPAY11 = [ 0 ];
GROUPAZ11 = [ 0 ];
GROUPANOM12 = [ 20 ];
GROUPAX12 = [ 0 ];
GROUPAY12 = [ 0 ];
GROUPAZ12 = [ 0 ];
GROUPANOM13 = [ 20 ];
GROUPAX13 = [ 0 ];
GROUPAY13 = [ 0 ];
GROUPAZ13 = [ 0 ];
GROUPANOM14 = [ 20 ];
GROUPAX14 = [ 0 ];
GROUPAY14 = [ 0 ];
GROUPAZ14 = [ 0 ];
GROUPANOM15 = [ 20 ];
GROUPAX15 = [ 0 ];
GROUPAY15 = [ 0 ];
GROUPAZ15 = [ 0 ];
GROUPBNOM1 = [ 1,1,8 ];
GROUPBX1 = [ -.597112,.895669,0 ];
GROUPBY1 = [ 0,0,0 ];
GROUPBZ1 = [ 2.638889,2.216667,1.900000 ];
GROUPBNOM2 = [ 1,1,8 ];
GROUPBX2 = [ -.597112,.895669,0 ];
GROUPBY2 = [ 0,0,0 ];
GROUPBZ2 = [ 2.738889,2.316667,2.000000 ];
GROUPBNOM3 = [ 1,1,8 ];
GROUPBX3 = [ -.597112,.895669,0 ];
GROUPBY3 = [ 0,0,0 ];
GROUPBZ3 = [ 2.838889,2.416667,2.100000 ];
GROUPBNOM4 = [ 1,1,8 ];
GROUPBX4 = [ -.597112,.895669,0 ];
GROUPBY4 = [ 0,0,0 ];
GROUPBZ4 = [ 2.938889,2.516667,2.200000 ];
GROUPBNOM5 = [ 1,1,8 ];
GROUPBX5 = [ -.597112,.895669,0 ];
GROUPBY5 = [ 0,0,0 ];
GROUPBZ5 = [ 3.038889,2.616667,2.300000 ];
GROUPBNOM6 = [ 1,1,8 ];
GROUPBX6 = [ -.597112,.895669,0 ];
GROUPBY6 = [ 0,0,0 ];
GROUPBZ6 = [ 3.138889,2.716667,2.400000 ];
GROUPBNOM7 = [ 1,1,8 ];
GROUPBX7 = [ -.597112,.895669,0 ];
GROUPBY7 = [ 0,0,0 ];
GROUPBZ7 = [ 3.238889,2.816667,2.500000 ];
GROUPBNOM8 = [ 1,1,8 ];
GROUPBX8 = [ -.597112,.895669,0 ];
GROUPBY8 = [ 0,0,0 ];
GROUPBZ8 = [ 3.338889,2.916667,2.600000 ];
GROUPBNOM9 = [ 1,1,8 ];
GROUPBX9 = [ -.597112,.895669,0 ];
GROUPBY9 = [ 0,0,0 ];
GROUPBZ9 = [ 3.438889,3.016667,2.700000 ];
GROUPBNOM10 = [ 1,1,8 ];
GROUPBX10 = [ -.597112,.895669,0 ];
GROUPBY10 = [ 0,0,0 ];
GROUPBZ10 = [ 3.538889,3.116667,2.800000 ];
GROUPBNOM11 = [ 1,1,8 ];
GROUPBX11 = [ -.597112,.895669,0 ];
GROUPBY11 = [ 0,0,0 ];
GROUPBZ11 = [ 3.638889,3.216667,2.900000 ];
GROUPBNOM12 = [ 1,1,8 ];
GROUPBX12 = [ -.597112,.895669,0 ];
GROUPBY12 = [ 0,0,0 ];
GROUPBZ12 = [ 4.238889,3.816667,3.500000 ];
GROUPBNOM13 = [ 1,1,8 ];
GROUPBX13 = [ -.597112,.895669,0 ];
GROUPBY13 = [ 0,0,0 ];
GROUPBZ13 = [ 4.738889,4.316667,4.000000 ];
GROUPBNOM14 = [ 1,1,8 ];
GROUPBX14 = [ -.597112,.895669,0 ];
GROUPBY14 = [ 0,0,0 ];
GROUPBZ14 = [ 5.238889,4.816667,4.500000 ];
GROUPBNOM15 = [ 1,1,8 ];
GROUPBX15 = [ -.597112,.895669,0 ];
GROUPBY15 = [ 0,0,0 ];
GROUPBZ15 = [ 5.738889,5.316667,5.000000 ];
DELTA = [ 
100.03,65.9,43.63,28.9,19.07,12.48,8.05,5.07,3.08,1.76,0.9,-0.34,-0.2,-0.09,-0.04
 ];
RTEMP = 
[sqrt((GROUPAX1(1)-GROUPBX1(2))^2+(GROUPAY1(1)-GROUPBY1(2))^2+(GROUPAZ1(1)-GROUPBZ1(2))^2),sqrt((GROUPAX2(1)-GROUPBX2(2))^2+(GROUPAY2(1)-GROUPBY2(2))^2+(GROUPAZ2(1)-GROUPBZ2(2))^2),sqrt((GROUPAX3(1)-GROUPBX3(2))^2+(GROUPAY3(1)-GROUPBY3(2))^2+(GROUPAZ3(1)-GROUPBZ3(2))^2),sqrt((GROUPAX4(1)-GROUPBX4(2))^2+(GROUPAY4(1)-GROUPBY4(2))^2+(GROUPAZ4(1)-GROUPBZ4(2))^2),sqrt((GROUPAX5(1)-GROUPBX5(2))^2+(GROUPAY5(1)-GROUPBY5(2))^2+(GROUPAZ5(1)-GROUPBZ5(2))^2),sqrt((GROUPAX6(1)-GROUPBX6(2))^2+(GROUPAY6(1)-GROUPBY6(2))^2+(GROUPAZ6(1)-GROUPBZ6(2))^2),sqrt((GROUPAX7(1)-GROUPBX7(2))^2+(GROUPAY7(1)-GROUPBY7(2))^2+(GROUPAZ7(1)-GROUPBZ7(2))^2),sqrt((GROUPAX8(1)-GROUPBX8(2))^2+(GROUPAY8(1)-GROUPBY8(2))^2+(GROUPAZ8(1)-GROUPBZ8(2))^2),sqrt((GROUPAX9(1)-GROUPBX9(2))^2+(GROUPAY9(1)-GROUPBY9(2))^2+(GROUPAZ9(1)-GROUPBZ9(2))^2),sqrt((GROUPAX10(1)-GROUPBX10(2))^2+(GROUPAY10(1)-GROUPBY10(2))^2+(GROUPAZ10(1)-GROUPBZ10(2))^2),sqrt((GROUPAX11(1)-GROUPBX11(2))^2+(GROUPAY11(1)-GROUPBY11(2))^2+(GROUPAZ11(1)-GROUPBZ11(2))^2),sqrt((GROUPAX12(1)-GROUPBX12(2))^2+(GROUPAY12(1)-GROUPBY12(2))^2+(GROUPAZ12(1)-GROUPBZ12(2))^2),sqrt((GROUPAX13(1)-GROUPBX13(2))^2+(GROUPAY13(1)-GROUPBY13(2))^2+(GROUPAZ13(1)-GROUPBZ13(2))^2),sqrt((GROUPAX14(1)-GROUPBX14(2))^2+(GROUPAY14(1)-GROUPBY14(2))^2+(GROUPAZ14(1)-GROUPBZ14(2))^2),sqrt((GROUPAX15(1)-GROUPBX15(2))^2+(GROUPAY15(1)-GROUPBY15(2))^2+(GROUPAZ15(1)-GROUPBZ15(2))^2);sqrt((GROUPAX1(1)-GROUPBX1(2))^2+(GROUPAY1(1)-GROUPBY1(2))^2+(GROUPAZ1(1)-GROUPBZ1(2))^2),sqrt((GROUPAX2(1)-GROUPBX2(2))^2+(GROUPAY2(1)-GROUPBY2(2))^2+(GROUPAZ2(1)-GROUPBZ2(2))^2),sqrt((GROUPAX3(1)-GROUPBX3(2))^2+(GROUPAY3(1)-GROUPBY3(2))^2+(GROUPAZ3(1)-GROUPBZ3(2))^2),sqrt((GROUPAX4(1)-GROUPBX4(2))^2+(GROUPAY4(1)-GROUPBY4(2))^2+(GROUPAZ4(1)-GROUPBZ4(2))^2),sqrt((GROUPAX5(1)-GROUPBX5(2))^2+(GROUPAY5(1)-GROUPBY5(2))^2+(GROUPAZ5(1)-GROUPBZ5(2))^2),sqrt((GROUPAX6(1)-GROUPBX6(2))^2+(GROUPAY6(1)-GROUPBY6(2))^2+(GROUPAZ6(1)-GROUPBZ6(2))^2),sqrt((GROUPAX7(1)-GROUPBX7(2))^2+(GROUPAY7(1)-GROUPBY7(2))^2+(GROUPAZ7(1)-GROUPBZ7(2))^2),sqrt((GROUPAX8(1)-GROUPBX8(2))^2+(GROUPAY8(1)-GROUPBY8(2))^2+(GROUPAZ8(1)-GROUPBZ8(2))^2),sqrt((GROUPAX9(1)-GROUPBX9(2))^2+(GROUPAY9(1)-GROUPBY9(2))^2+(GROUPAZ9(1)-GROUPBZ9(2))^2),sqrt((GROUPAX10(1)-GROUPBX10(2))^2+(GROUPAY10(1)-GROUPBY10(2))^2+(GROUPAZ10(1)-GROUPBZ10(2))^2),sqrt((GROUPAX11(1)-GROUPBX11(2))^2+(GROUPAY11(1)-GROUPBY11(2))^2+(GROUPAZ11(1)-GROUPBZ11(2))^2),sqrt((GROUPAX12(1)-GROUPBX12(2))^2+(GROUPAY12(1)-GROUPBY12(2))^2+(GROUPAZ12(1)-GROUPBZ12(2))^2),sqrt((GROUPAX13(1)-GROUPBX13(2))^2+(GROUPAY13(1)-GROUPBY13(2))^2+(GROUPAZ13(1)-GROUPBZ13(2))^2),sqrt((GROUPAX14(1)-GROUPBX14(2))^2+(GROUPAY14(1)-GROUPBY14(2))^2+(GROUPAZ14(1)-GROUPBZ14(2))^2),sqrt((GROUPAX15(1)-GROUPBX15(2))^2+(GROUPAY15(1)-GROUPBY15(2))^2+(GROUPAZ15(1)-GROUPBZ15(2))^2);sqrt((GROUPAX1(1)-GROUPBX1(3))^2+(GROUPAY1(1)-GROUPBY1(3))^2+(GROUPAZ1(1)-GROUPBZ1(3))^2),sqrt((GROUPAX2(1)-GROUPBX2(3))^2+(GROUPAY2(1)-GROUPBY2(3))^2+(GROUPAZ2(1)-GROUPBZ2(3))^2),sqrt((GROUPAX3(1)-GROUPBX3(3))^2+(GROUPAY3(1)-GROUPBY3(3))^2+(GROUPAZ3(1)-GROUPBZ3(3))^2),sqrt((GROUPAX4(1)-GROUPBX4(3))^2+(GROUPAY4(1)-GROUPBY4(3))^2+(GROUPAZ4(1)-GROUPBZ4(3))^2),sqrt((GROUPAX5(1)-GROUPBX5(3))^2+(GROUPAY5(1)-GROUPBY5(3))^2+(GROUPAZ5(1)-GROUPBZ5(3))^2),sqrt((GROUPAX6(1)-GROUPBX6(3))^2+(GROUPAY6(1)-GROUPBY6(3))^2+(GROUPAZ6(1)-GROUPBZ6(3))^2),sqrt((GROUPAX7(1)-GROUPBX7(3))^2+(GROUPAY7(1)-GROUPBY7(3))^2+(GROUPAZ7(1)-GROUPBZ7(3))^2),sqrt((GROUPAX8(1)-GROUPBX8(3))^2+(GROUPAY8(1)-GROUPBY8(3))^2+(GROUPAZ8(1)-GROUPBZ8(3))^2),sqrt((GROUPAX9(1)-GROUPBX9(3))^2+(GROUPAY9(1)-GROUPBY9(3))^2+(GROUPAZ9(1)-GROUPBZ9(3))^2),sqrt((GROUPAX10(1)-GROUPBX10(3))^2+(GROUPAY10(1)-GROUPBY10(3))^2+(GROUPAZ10(1)-GROUPBZ10(3))^2),sqrt((GROUPAX11(1)-GROUPBX11(3))^2+(GROUPAY11(1)-GROUPBY11(3))^2+(GROUPAZ11(1)-GROUPBZ11(3))^2),sqrt((GROUPAX12(1)-GROUPBX12(3))^2+(GROUPAY12(1)-GROUPBY12(3))^2+(GROUPAZ12(1)-GROUPBZ12(3))^2),sqrt((GROUPAX13(1)-GROUPBX13(3))^2+(GROUPAY13(1)-GROUPBY13(3))^2+(GROUPAZ13(1)-GROUPBZ13(3))^2),sqrt((GROUPAX14(1)-GROUPBX14(3))^2+(GROUPAY14(1)-GROUPBY14(3))^2+(GROUPAZ14(1)-GROUPBZ14(3))^2),sqrt((GROUPAX15(1)-GROUPBX15(3))^2+(GROUPAY15(1)-GROUPBY15(3))^2+(GROUPAZ15(1)-GROUPBZ15(3))^2)
 ];