src/math/SquareMatrix.hpp

/*
 * Copyright (C) 2000-2004  Object Oriented Parallel Simulation Engine (OOPSE) project
 * 
 * Contact: oopse@oopse.org
 * 
 * This program is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public License
 * as published by the Free Software Foundation; either version 2.1
 * of the License, or (at your option) any later version.
 * All we ask is that proper credit is given for our work, which includes
 * - but is not limited to - adding the above copyright notice to the beginning
 * of your source code files, and to any copyright notice that you may distribute
 * with programs based on this work.
 * 
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU Lesser General Public License for more details.
 * 
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
 *
 */

/**
 * @file SquareMatrix.hpp
 * @author Teng Lin
 * @date 10/11/2004
 * @version 1.0
 */
 #ifndef MATH_SQUAREMATRIX_HPP 
#define MATH_SQUAREMATRIX_HPP 

#include "math/RectMatrix.hpp"

namespace oopse {

    /**
     * @class SquareMatrix SquareMatrix.hpp "math/SquareMatrix.hpp"
     * @brief A square matrix class
     * @template Real the element type
     * @template Dim the dimension of the square matrix
     */
    template<typename Real, int Dim>
    class SquareMatrix : public RectMatrix<Real, Dim, Dim> {
        public:
            typedef Real ElemType;
            typedef Real* ElemPoinerType;

            /** default constructor */
            SquareMatrix() {
                for (unsigned int i = 0; i < Dim; i++)
                    for (unsigned int j = 0; j < Dim; j++)
                        data_[i][j] = 0.0;
             }

            /** copy constructor */
            SquareMatrix(const RectMatrix<Real, Dim, Dim>& m) : RectMatrix<Real, Dim, Dim>(m) {
            }
            
            /** copy assignment operator */
            SquareMatrix<Real, Dim>& operator =(const RectMatrix<Real, Dim, Dim>& m) {
                RectMatrix<Real, Dim, Dim>::operator=(m);
                return *this;
            }
                                   
            /** Retunrs  an identity matrix*/

           static SquareMatrix<Real, Dim> identity() {
                SquareMatrix<Real, Dim> m;
                
                for (unsigned int i = 0; i < Dim; i++) 
                    for (unsigned int j = 0; j < Dim; j++) 
                        if (i == j)
                            m(i, j) = 1.0;
                        else
                            m(i, j) = 0.0;

                return m;
            }

            /** 
             * Retunrs  the inversion of this matrix. 
             * @todo need implementation
             */
             SquareMatrix<Real, Dim>  inverse() {
                 SquareMatrix<Real, Dim> result;

                 return result;
            }        

            /**
             * Returns the determinant of this matrix.
             * @todo need implementation
             */
            Real determinant() const {
                Real det;
                return det;
            }

            /** Returns the trace of this matrix. */
            Real trace() const {
               Real tmp = 0;
               
                for (unsigned int i = 0; i < Dim ; i++)
                    tmp += data_[i][i];

                return tmp;
            }

            /** Tests if this matrix is symmetrix. */            
            bool isSymmetric() const {
                for (unsigned int i = 0; i < Dim - 1; i++)
                    for (unsigned int j = i; j < Dim; j++)
                        if (fabs(data_[i][j] - data_[j][i]) > oopse::epsilon) 
                            return false;
                        
                return true;
            }

            /** Tests if this matrix is orthogonal. */            
            bool isOrthogonal() {
                SquareMatrix<Real, Dim> tmp;

                tmp = *this * transpose();

                return tmp.isDiagonal();
            }

            /** Tests if this matrix is diagonal. */
            bool isDiagonal() const {
                for (unsigned int i = 0; i < Dim ; i++)
                    for (unsigned int j = 0; j < Dim; j++)
                        if (i !=j && fabs(data_[i][j]) > oopse::epsilon) 
                            return false;
                        
                return true;
            }

            /** Tests if this matrix is the unit matrix. */
            bool isUnitMatrix() const {
                if (!isDiagonal())
                    return false;
                
                for (unsigned int i = 0; i < Dim ; i++)
                    if (fabs(data_[i][i] - 1) > oopse::epsilon)
                        return false;
                    
                return true;
            }         

            /** @todo need implementation */
            void diagonalize() {
                //jacobi(m, eigenValues, ortMat);
            }

            /**
             * Jacobi iteration routines for computing eigenvalues/eigenvectors of 
             * real symmetric matrix
             *
             * @return true if success, otherwise return false
             * @param a symmetric matrix whose eigenvectors are to be computed. On return, the matrix is
             *     overwritten
             * @param w will contain the eigenvalues of the matrix On return of this function
             * @param v the columns of this matrix will contain the eigenvectors. The eigenvectors are 
             *    normalized and mutually orthogonal. 
             */
           
            static int jacobi(SquareMatrix<Real, Dim>& a, Vector<Real, Dim>& d, 
                                  SquareMatrix<Real, Dim>& v);
    };//end SquareMatrix


/*=========================================================================

  Program:   Visualization Toolkit
  Module:    $RCSfile: SquareMatrix.hpp,v $

  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
     PURPOSE.  See the above copyright notice for more information.

=========================================================================*/

#define VTK_ROTATE(a,i,j,k,l) g=a(i, j);h=a(k, l);a(i, j)=g-s*(h+g*tau);\
        a(k, l)=h+s*(g-h*tau)

#define VTK_MAX_ROTATIONS 20

    // Jacobi iteration for the solution of eigenvectors/eigenvalues of a nxn
    // real symmetric matrix. Square nxn matrix a; size of matrix in n;
    // output eigenvalues in w; and output eigenvectors in v. Resulting
    // eigenvalues/vectors are sorted in decreasing order; eigenvectors are
    // normalized.
    template<typename Real, int Dim>
    int SquareMatrix<Real, Dim>::jacobi(SquareMatrix<Real, Dim>& a, Vector<Real, Dim>& w, 
                                  SquareMatrix<Real, Dim>& v) {
      const int n = Dim;  
      int i, j, k, iq, ip, numPos;
      Real tresh, theta, tau, t, sm, s, h, g, c, tmp;
      Real bspace[4], zspace[4];
      Real *b = bspace;
      Real *z = zspace;

      // only allocate memory if the matrix is large
      if (n > 4)
        {
        b = new Real[n];
        z = new Real[n]; 
        }

      // initialize
      for (ip=0; ip<n; ip++) 
        {
        for (iq=0; iq<n; iq++)
          {
          v(ip, iq) = 0.0;
          }
        v(ip, ip) = 1.0;
        }
      for (ip=0; ip<n; ip++) 
        {
        b[ip] = w[ip] = a(ip, ip);
        z[ip] = 0.0;
        }

      // begin rotation sequence
      for (i=0; i<VTK_MAX_ROTATIONS; i++) 
        {
        sm = 0.0;
        for (ip=0; ip<n-1; ip++) 
          {
          for (iq=ip+1; iq<n; iq++)
            {
            sm += fabs(a(ip, iq));
            }
          }
        if (sm == 0.0)
          {
          break;
          }

        if (i < 3)                                // first 3 sweeps
          {
          tresh = 0.2*sm/(n*n);
          }
        else
          {
          tresh = 0.0;
          }

        for (ip=0; ip<n-1; ip++) 
          {
          for (iq=ip+1; iq<n; iq++) 
            {
            g = 100.0*fabs(a(ip, iq));

            // after 4 sweeps
            if (i > 3 && (fabs(w[ip])+g) == fabs(w[ip])
            && (fabs(w[iq])+g) == fabs(w[iq]))
              {
              a(ip, iq) = 0.0;
              }
            else if (fabs(a(ip, iq)) > tresh) 
              {
              h = w[iq] - w[ip];
              if ( (fabs(h)+g) == fabs(h))
                {
                t = (a(ip, iq)) / h;
                }
              else 
                {
                theta = 0.5*h / (a(ip, iq));
                t = 1.0 / (fabs(theta)+sqrt(1.0+theta*theta));
                if (theta < 0.0)
                  {
                  t = -t;
                  }
                }
              c = 1.0 / sqrt(1+t*t);
              s = t*c;
              tau = s/(1.0+c);
              h = t*a(ip, iq);
              z[ip] -= h;
              z[iq] += h;
              w[ip] -= h;
              w[iq] += h;
              a(ip, iq)=0.0;

              // ip already shifted left by 1 unit
              for (j = 0;j <= ip-1;j++) 
                {
                VTK_ROTATE(a,j,ip,j,iq);
                }
              // ip and iq already shifted left by 1 unit
              for (j = ip+1;j <= iq-1;j++) 
                {
                VTK_ROTATE(a,ip,j,j,iq);
                }
              // iq already shifted left by 1 unit
              for (j=iq+1; j<n; j++) 
                {
                VTK_ROTATE(a,ip,j,iq,j);
                }
              for (j=0; j<n; j++) 
                {
                VTK_ROTATE(v,j,ip,j,iq);
                }
              }
            }
          }

        for (ip=0; ip<n; ip++) 
          {
          b[ip] += z[ip];
          w[ip] = b[ip];
          z[ip] = 0.0;
          }
        }

      //// this is NEVER called
      if ( i >= VTK_MAX_ROTATIONS )
        {
           std::cout << "vtkMath::Jacobi: Error extracting eigenfunctions" << std::endl;
           return 0;
        }

      // sort eigenfunctions                 these changes do not affect accuracy 
      for (j=0; j<n-1; j++)                  // boundary incorrect
        {
        k = j;
        tmp = w[k];
        for (i=j+1; i<n; i++)                // boundary incorrect, shifted already
          {
          if (w[i] >= tmp)                   // why exchage if same?
            {
            k = i;
            tmp = w[k];
            }
          }
        if (k != j) 
          {
          w[k] = w[j];
          w[j] = tmp;
          for (i=0; i<n; i++) 
            {
            tmp = v(i, j);
            v(i, j) = v(i, k);
            v(i, k) = tmp;
            }
          }
        }
      // insure eigenvector consistency (i.e., Jacobi can compute vectors that
      // are negative of one another (.707,.707,0) and (-.707,-.707,0). This can
      // reek havoc in hyperstreamline/other stuff. We will select the most
      // positive eigenvector.
      int ceil_half_n = (n >> 1) + (n & 1);
      for (j=0; j<n; j++)
        {
        for (numPos=0, i=0; i<n; i++)
          {
          if ( v(i, j) >= 0.0 )
            {
            numPos++;
            }
          }
    //    if ( numPos < ceil(double(n)/double(2.0)) )
        if ( numPos < ceil_half_n)
          {
          for(i=0; i<n; i++)
            {
            v(i, j) *= -1.0;
            }
          }
        }

      if (n > 4)
        {
        delete [] b;
        delete [] z;
        }
      return 1;
    }


}
#endif //MATH_SQUAREMATRIX_HPP 

Revision:	137
Committed:	Thu Oct 21 21:31:39 2004 UTC (21 years ago) by tim
File size:	12083 byte(s)
Log Message:	Snapshot and SnapshotManager in progress
#	User	Rev	Content
1	tim	70	/*
2			* Copyright (C) 2000-2004 Object Oriented Parallel Simulation Engine (OOPSE) project
3			*
4			* Contact: oopse@oopse.org
5			*
6			* This program is free software; you can redistribute it and/or
7			* modify it under the terms of the GNU Lesser General Public License
8			* as published by the Free Software Foundation; either version 2.1
9			* of the License, or (at your option) any later version.
10			* All we ask is that proper credit is given for our work, which includes
11			* - but is not limited to - adding the above copyright notice to the beginning
12			* of your source code files, and to any copyright notice that you may distribute
13			* with programs based on this work.
14			*
15			* This program is distributed in the hope that it will be useful,
16			* but WITHOUT ANY WARRANTY; without even the implied warranty of
17			* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18			* GNU Lesser General Public License for more details.
19			*
20			* You should have received a copy of the GNU Lesser General Public License
21			* along with this program; if not, write to the Free Software
22			* Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
23			*
24			*/
25
26			/**
27			* @file SquareMatrix.hpp
28			* @author Teng Lin
29			* @date 10/11/2004
30			* @version 1.0
31			*/
32	tim	123	#ifndef MATH_SQUAREMATRIX_HPP
33	tim	70	#define MATH_SQUAREMATRIX_HPP
34
35	tim	74	#include "math/RectMatrix.hpp"
36	tim	70
37			namespace oopse {
38
39			/**
40			* @class SquareMatrix SquareMatrix.hpp "math/SquareMatrix.hpp"
41			* @brief A square matrix class
42			* @template Real the element type
43			* @template Dim the dimension of the square matrix
44			*/
45			template<typename Real, int Dim>
46	tim	74	class SquareMatrix : public RectMatrix<Real, Dim, Dim> {
47	tim	70	public:
48	tim	137	typedef Real ElemType;
49			typedef Real* ElemPoinerType;
50	tim	70
51	tim	137	/** default constructor */
52			SquareMatrix() {
53			for (unsigned int i = 0; i < Dim; i++)
54			for (unsigned int j = 0; j < Dim; j++)
55			data_[i][j] = 0.0;
56			}
57	tim	70
58	tim	137	/** copy constructor */
59			SquareMatrix(const RectMatrix<Real, Dim, Dim>& m) : RectMatrix<Real, Dim, Dim>(m) {
60			}
61	tim	70
62	tim	137	/** copy assignment operator */
63			SquareMatrix<Real, Dim>& operator =(const RectMatrix<Real, Dim, Dim>& m) {
64			RectMatrix<Real, Dim, Dim>::operator=(m);
65			return *this;
66			}
67
68			/** Retunrs an identity matrix*/
69	tim	74
70	tim	137	static SquareMatrix<Real, Dim> identity() {
71			SquareMatrix<Real, Dim> m;
72
73			for (unsigned int i = 0; i < Dim; i++)
74			for (unsigned int j = 0; j < Dim; j++)
75			if (i == j)
76			m(i, j) = 1.0;
77			else
78			m(i, j) = 0.0;
79	tim	70
80	tim	137	return m;
81			}
82	tim	74
83	tim	137	/**
84			* Retunrs the inversion of this matrix.
85			* @todo need implementation
86			*/
87			SquareMatrix<Real, Dim> inverse() {
88			SquareMatrix<Real, Dim> result;
89	tim	70
90	tim	137	return result;
91			}
92	tim	70
93	tim	137	/**
94			* Returns the determinant of this matrix.
95			* @todo need implementation
96			*/
97			Real determinant() const {
98			Real det;
99			return det;
100			}
101	tim	70
102	tim	137	/** Returns the trace of this matrix. */
103			Real trace() const {
104			Real tmp = 0;
105
106			for (unsigned int i = 0; i < Dim ; i++)
107			tmp += data_[i][i];
108	tim	70
109	tim	137	return tmp;
110			}
111	tim	70
112	tim	137	/** Tests if this matrix is symmetrix. */
113			bool isSymmetric() const {
114			for (unsigned int i = 0; i < Dim - 1; i++)
115			for (unsigned int j = i; j < Dim; j++)
116			if (fabs(data_[i][j] - data_[j][i]) > oopse::epsilon)
117			return false;
118
119			return true;
120			}
121	tim	70
122	tim	137	/** Tests if this matrix is orthogonal. */
123			bool isOrthogonal() {
124			SquareMatrix<Real, Dim> tmp;
125	tim	70
126	tim	137	tmp = this transpose();
127	tim	70
128	tim	137	return tmp.isDiagonal();
129			}
130	tim	70
131	tim	137	/** Tests if this matrix is diagonal. */
132			bool isDiagonal() const {
133			for (unsigned int i = 0; i < Dim ; i++)
134			for (unsigned int j = 0; j < Dim; j++)
135			if (i !=j && fabs(data_[i][j]) > oopse::epsilon)
136			return false;
137
138			return true;
139			}
140
141			/** Tests if this matrix is the unit matrix. */
142			bool isUnitMatrix() const {
143			if (!isDiagonal())
144	tim	70	return false;
145
146	tim	137	for (unsigned int i = 0; i < Dim ; i++)
147			if (fabs(data_[i][i] - 1) > oopse::epsilon)
148			return false;
149
150			return true;
151			}
152	tim	70
153	tim	137	/** @todo need implementation */
154			void diagonalize() {
155			//jacobi(m, eigenValues, ortMat);
156			}
157	tim	76
158	tim	137	/**
159			* Jacobi iteration routines for computing eigenvalues/eigenvectors of
160			* real symmetric matrix
161			*
162			* @return true if success, otherwise return false
163			* @param a symmetric matrix whose eigenvectors are to be computed. On return, the matrix is
164			* overwritten
165			* @param w will contain the eigenvalues of the matrix On return of this function
166			* @param v the columns of this matrix will contain the eigenvectors. The eigenvectors are
167			* normalized and mutually orthogonal.
168			*/
169
170			static int jacobi(SquareMatrix<Real, Dim>& a, Vector<Real, Dim>& d,
171			SquareMatrix<Real, Dim>& v);
172	tim	70	};//end SquareMatrix
173
174	tim	76
175	tim	123	/*=========================================================================
176	tim	76
177	tim	123	Program: Visualization Toolkit
178			Module: $RCSfile: SquareMatrix.hpp,v $
179	tim	76
180	tim	123	Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
181			All rights reserved.
182			See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
183
184			This software is distributed WITHOUT ANY WARRANTY; without even
185			the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
186			PURPOSE. See the above copyright notice for more information.
187
188			=========================================================================*/
189
190			#define VTK_ROTATE(a,i,j,k,l) g=a(i, j);h=a(k, l);a(i, j)=g-s(h+gtau);\
191			a(k, l)=h+s(g-htau)
192
193			#define VTK_MAX_ROTATIONS 20
194
195			// Jacobi iteration for the solution of eigenvectors/eigenvalues of a nxn
196			// real symmetric matrix. Square nxn matrix a; size of matrix in n;
197			// output eigenvalues in w; and output eigenvectors in v. Resulting
198			// eigenvalues/vectors are sorted in decreasing order; eigenvectors are
199			// normalized.
200			template<typename Real, int Dim>
201			int SquareMatrix<Real, Dim>::jacobi(SquareMatrix<Real, Dim>& a, Vector<Real, Dim>& w,
202			SquareMatrix<Real, Dim>& v) {
203			const int n = Dim;
204			int i, j, k, iq, ip, numPos;
205			Real tresh, theta, tau, t, sm, s, h, g, c, tmp;
206			Real bspace[4], zspace[4];
207			Real *b = bspace;
208			Real *z = zspace;
209
210			// only allocate memory if the matrix is large
211			if (n > 4)
212			{
213			b = new Real[n];
214			z = new Real[n];
215			}
216
217			// initialize
218			for (ip=0; ip<n; ip++)
219			{
220			for (iq=0; iq<n; iq++)
221			{
222			v(ip, iq) = 0.0;
223			}
224	tim	93	v(ip, ip) = 1.0;
225	tim	123	}
226			for (ip=0; ip<n; ip++)
227			{
228			b[ip] = w[ip] = a(ip, ip);
229			z[ip] = 0.0;
230			}
231	tim	76
232	tim	123	// begin rotation sequence
233			for (i=0; i<VTK_MAX_ROTATIONS; i++)
234			{
235	tim	93	sm = 0.0;
236	tim	123	for (ip=0; ip<n-1; ip++)
237			{
238			for (iq=ip+1; iq<n; iq++)
239			{
240			sm += fabs(a(ip, iq));
241			}
242			}
243	tim	93	if (sm == 0.0)
244	tim	123	{
245			break;
246			}
247	tim	76
248	tim	123	if (i < 3) // first 3 sweeps
249			{
250			tresh = 0.2sm/(nn);
251			}
252	tim	93	else
253	tim	123	{
254			tresh = 0.0;
255			}
256	tim	76
257	tim	123	for (ip=0; ip<n-1; ip++)
258			{
259			for (iq=ip+1; iq<n; iq++)
260			{
261			g = 100.0*fabs(a(ip, iq));
262	tim	76
263	tim	123	// after 4 sweeps
264			if (i > 3 && (fabs(w[ip])+g) == fabs(w[ip])
265			&& (fabs(w[iq])+g) == fabs(w[iq]))
266			{
267			a(ip, iq) = 0.0;
268			}
269			else if (fabs(a(ip, iq)) > tresh)
270			{
271			h = w[iq] - w[ip];
272			if ( (fabs(h)+g) == fabs(h))
273			{
274			t = (a(ip, iq)) / h;
275			}
276			else
277			{
278			theta = 0.5*h / (a(ip, iq));
279			t = 1.0 / (fabs(theta)+sqrt(1.0+theta*theta));
280			if (theta < 0.0)
281			{
282			t = -t;
283			}
284			}
285			c = 1.0 / sqrt(1+t*t);
286			s = t*c;
287			tau = s/(1.0+c);
288			h = t*a(ip, iq);
289			z[ip] -= h;
290			z[iq] += h;
291			w[ip] -= h;
292			w[iq] += h;
293			a(ip, iq)=0.0;
294	tim	76
295	tim	123	// ip already shifted left by 1 unit
296			for (j = 0;j <= ip-1;j++)
297			{
298			VTK_ROTATE(a,j,ip,j,iq);
299	tim	93	}
300	tim	123	// ip and iq already shifted left by 1 unit
301			for (j = ip+1;j <= iq-1;j++)
302			{
303			VTK_ROTATE(a,ip,j,j,iq);
304			}
305			// iq already shifted left by 1 unit
306			for (j=iq+1; j<n; j++)
307			{
308			VTK_ROTATE(a,ip,j,iq,j);
309			}
310			for (j=0; j<n; j++)
311			{
312			VTK_ROTATE(v,j,ip,j,iq);
313			}
314			}
315	tim	93	}
316	tim	123	}
317	tim	93
318	tim	123	for (ip=0; ip<n; ip++)
319			{
320			b[ip] += z[ip];
321			w[ip] = b[ip];
322			z[ip] = 0.0;
323			}
324	tim	93	}
325
326	tim	123	//// this is NEVER called
327			if ( i >= VTK_MAX_ROTATIONS )
328			{
329			std::cout << "vtkMath::Jacobi: Error extracting eigenfunctions" << std::endl;
330			return 0;
331			}
332	tim	76
333	tim	123	// sort eigenfunctions these changes do not affect accuracy
334			for (j=0; j<n-1; j++) // boundary incorrect
335			{
336	tim	93	k = j;
337	tim	123	tmp = w[k];
338			for (i=j+1; i<n; i++) // boundary incorrect, shifted already
339			{
340			if (w[i] >= tmp) // why exchage if same?
341			{
342	tim	93	k = i;
343	tim	123	tmp = w[k];
344	tim	93	}
345	tim	123	}
346			if (k != j)
347			{
348			w[k] = w[j];
349			w[j] = tmp;
350			for (i=0; i<n; i++)
351			{
352			tmp = v(i, j);
353			v(i, j) = v(i, k);
354			v(i, k) = tmp;
355			}
356			}
357	tim	93	}
358	tim	123	// insure eigenvector consistency (i.e., Jacobi can compute vectors that
359			// are negative of one another (.707,.707,0) and (-.707,-.707,0). This can
360			// reek havoc in hyperstreamline/other stuff. We will select the most
361			// positive eigenvector.
362			int ceil_half_n = (n >> 1) + (n & 1);
363			for (j=0; j<n; j++)
364			{
365			for (numPos=0, i=0; i<n; i++)
366			{
367			if ( v(i, j) >= 0.0 )
368			{
369			numPos++;
370	tim	93	}
371	tim	123	}
372			// if ( numPos < ceil(double(n)/double(2.0)) )
373			if ( numPos < ceil_half_n)
374			{
375			for(i=0; i<n; i++)
376			{
377			v(i, j) *= -1.0;
378			}
379			}
380	tim	93	}
381	tim	76
382	tim	123	if (n > 4)
383			{
384			delete [] b;
385			delete [] z;
386			}
387			return 1;
388	tim	76	}
389
390
391			}
392	tim	123	#endif //MATH_SQUAREMATRIX_HPP
393	tim	76