Метод Гаусса - ошибка Floating point - C (СИ)

Формулировка задачи:

Borland C версия 3.1 Решение системы методом Гаусса. Выскакивает ошибка Floating point:square root of negative number во время исполнения ф-ии stage2. Я убилась,но мне не отладить. Помогите,пожалуйста=))

#include <stdio.h>
#include <conio.h>
#include <math.h>
#define size 10
 
void input(int   n,int m,float A[size][size],float B[]);
void stage1(float A[size][size],float L[size][size],float U[size][size],int n);
void stage2(float X[size],float B[size],float L[size][size],int n);
void stage3(float X[size],float B[size],float U[size][size],int n);
void result (int n,float Z[]);
 
void main()
{
 int n,m,i,j,k;
 float A[size][size],B[size],X[size],L[size][size],U[size][size];
 puts("Input:");
 printf("Give the size of the matrix n_m:\n");
 scanf("%i%i",&n,&m);
 input(n,m,A,B);
 stage1(A,L,U,n); 
 stage2(X,B,L,n);
 stag[CPP][/CPP]e3(X,B,U,n);
 result(n,X);
   
}
 
 void input(int n,int m,float A[size][size],float B[])
 {
  int i,j; 
  puts("Enter matrix A[][]:");
  for(i=0; i<n; i++)
      { for (j=0; j<m; j++)
             scanf("%f",&A[i][j]);
      }
  puts("Enter vector B[]:");
  for (i=0;i<n;i++)    
       scanf("%f",&B[i]);      
 
 }

 void result (int n,float Z[])
{
 int i;
 puts("Vector X:");
 for (i=0;i<n;i++)
     printf("%5.2f\n",Z[i]);
}
 void stage1(float A[size][size],float L[size][size],float U[size][size],int n)
 { int i,j,k;
    for(i=0; i<n; i++)
      { for (j=0; j<n; j++)
        L[i][j]=0;
        U[i][j]=0;
       }
   for(i=0; i<n; i++)
      { for (j=0; j<n; j++)
         
         if(i>j) L[i][j]=A[i][j];
        else U[i][j]=A[i][j];
       }
   for (k=0;k<n-1;k++)
    { for(i=k+1; i<n; i++)  
      A[i][k]=A[i][k]/A[k][k];
          for (j=k+1; j<n; j++)
                  A[i][j]=A[i][j]-A[i][k]*A[k][j];
 
         }
 
 }
 void stage2(float X[size],float B[size],float L[size][size],int n)
 { int i,j;
 
   for(i=0;i<n;i++)
     { 
      X[i]=B[i];
        
    for (j=0;j<i-1;j++)    
        {
          X[i]=X[i]-L[i][j]*X[j];
              X[i]=X[i]/L[i][i];             
             }
     }
 }
 
  void stage3(float X[size],float B[size],float U[size][size],int n)
  { int i,j;

    for(i=n-1; i>=0; i--)
    { X[i]=B[i];
     for (j=i+1;j<n;j++)
         {X[i]=X[i] - U[i][j]*X[j];
          X[i]=X[i]/U[i][i];
              }
         }
  }

Решение задачи: «Метод Гаусса - ошибка Floating point»

textual
                            Листинг программы 
                        

                            #include "stdafx.h"
#include "iostream"
#include "cmath"
#include "conio.h"
 
using namespace std;#include "stdafx.h"
#include "iostream"
#include "cmath"
#include "conio.h"
 
using namespace std;
 
void OutMatrix(double **a, int n)
{
    for(int i = 0; i<n; i++)
    {
        for(int j = 0; j<n; j++)
            cout << a[i][j] << "    ";
        cout << endl;
    }
}
 
void OutVector(double *x, int n)
{
    for(int i = 0; i<n; i++)
            cout << x[i] <<"    ";
    cout << endl;
}
 
 
void Gauss(double **a, double *b, int n, double *x)
{
    // Прямой ход метода 
    double M;
    for(int i = 0; i<n-1; i++)
    {
        for(int j = i; j<n-1; j++)
        {
            M = -a[j+1][i]/a[i][i];
            for(int k = i; k<n; k++)
            {
                a[j+1][k] += a[i][k]*M; 
                if(a[j+1][k] < 0.000001 & a[j+1][k]> -0.000001)
                    a[j+1][k] = 0;
            }
            b[j+1] += b[i]*M;
        }
    }
    // Обратный ход 
    int N = n-1;
    for(int i = 0; i<=N; i++)
        for(int j = N; j>=N-i; j--)
        {
            if( j == N)
                x[N-i] = b[N-i];
            if(j>(N-i) & j<N)
                x[N-i] -= x[j];
            if(j == N-i)
                x[N-i] /= a[j][j];
        }
}
 
void MatrMultVect(double **a1, double *x, int n, double *b1)
{
    double b;
    for(int i = 0; i<n; i++)
    {
        b = 0;
        for(int j = 0; j<n; j++)
            b += a1[i][j]*x[j];
        b1[i] = b;
        
    }
}
 
void OutMatrix(double **a, int n)
{
    for(int i = 0; i<n; i++)
    {
        for(int j = 0; j<n; j++)
            cout << a[i][j] << "    ";
        cout << endl;
    }
}
 
void OutVector(double *x, int n)
{
    for(int i = 0; i<n; i++)
            cout << x[i] <<"    ";
    cout << endl;
}
 
 
void Gauss(double **a, double *b, int n, double *x)
{
    // Прямой ход метода 
    double M;
    for(int i = 0; i<n-1; i++)
    {
        for(int j = i; j<n-1; j++)
        {
            M = -a[j+1][i]/a[i][i];
            for(int k = i; k<n; k++)
            {
                a[j+1][k] += a[i][k]*M; 
                if(a[j+1][k] < 0.000001 & a[j+1][k]> -0.000001)
                    a[j+1][k] = 0;
            }
            b[j+1] += b[i]*M;
        }
    }
    // Обратный ход 
    int N = n-1;
    for(int i = 0; i<=N; i++)
        for(int j = N; j>=N-i; j--)
        {
            if( j == N)
                x[N-i] = b[N-i];
            if(j>(N-i) & j<N)
                x[N-i] -= x[j];
            if(j == N-i)
                x[N-i] /= a[j][j];
        }
}
 
void MatrMultVect(double **a1, double *x, int n, double *b1)
{
    double b;
    for(int i = 0; i<n; i++)
    {
        b = 0;
        for(int j = 0; j<n; j++)
            b += a1[i][j]*x[j];
        b1[i] = b;
        
    }
}
                        