/* NUMclasses.c
 *
 * Copyright (C) 1997-2002 David Weenink
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 *
 * 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
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/*
 djmw 20010719
 djmw 20020408 GPL + cosmetic changes
*/

#include "NUMclasses.h"
#include "NUMlapack.h"
#include "NUMmachar.h"
#include "Collection.h"

#include "oo_DESTROY.h"
#include "NUMclasses_def.h"
#include "oo_COPY.h"
#include "NUMclasses_def.h"
#include "oo_EQUAL.h"
#include "NUMclasses_def.h"
#include "oo_WRITE_ASCII.h"
#include "NUMclasses_def.h"
#include "oo_WRITE_BINARY.h"
#include "NUMclasses_def.h"
#include "oo_READ_ASCII.h"
#include "NUMclasses_def.h"
#include "oo_READ_BINARY.h"
#include "NUMclasses_def.h"
#include "oo_DESCRIPTION.h"
#include "NUMclasses_def.h"

extern machar_Table NUMfpp;

static void classSVD_info (I)
{
	iam (SVD);
	Melder_information ("Number of rows = %ld\n"
		"Number of columns = %ld\n",
		my numberOfRows, my numberOfColumns);
}

class_methods (SVD, Data)
	class_method_local (SVD, destroy)
	class_method_local (SVD, equal)
	class_method_local (SVD, copy)
	class_method_local (SVD, readAscii)
	class_method_local (SVD, readBinary)
	class_method_local (SVD, writeAscii)
	class_method_local (SVD, writeBinary)
	class_method_local (SVD, description)
	class_method_local (SVD, info)
class_methods_end

SVD SVD_create (long numberOfRows, long numberOfColumns)
{
	SVD me = new (SVD);
	if (! me) return NULL;
	my numberOfRows = numberOfRows;
	my numberOfColumns = numberOfColumns;
	if (! NUMfpp) NUMmachar ();
	my tolerance = NUMfpp -> eps_d * MAX (numberOfRows, numberOfColumns);
	if (((my u = NUMdmatrix (1, numberOfRows, 1, numberOfColumns)) == NULL) ||
		((my v = NUMdmatrix (1, numberOfColumns, 1, numberOfColumns)) == NULL) ||
		((my d = NUMdvector (1, numberOfColumns)) == NULL)) forget (me);
	return me;
}

SVD SVD_create_d (double **m, long numberOfRows, long numberOfColumns)
{
	SVD me = SVD_create (numberOfRows, numberOfColumns);
	
	if ((me == NULL) || ! SVD_svd_d (me, m)) forget (me);
	return me;
}

SVD SVD_create_f (float **m, long numberOfRows, long numberOfColumns)
{
	SVD me = SVD_create (numberOfRows, numberOfColumns);
	 
	if ((me == NULL) || ! SVD_svd_f (me, m)) forget (me);
	return me;
}

int SVD_svd_d (SVD me, double **m)
{
	NUMdmatrix_copyElements (m, my u, 1, my numberOfRows, 1, 
		my numberOfColumns);
	return NUMsvdcmp (my u, my numberOfRows, my numberOfColumns, my d, my v);
}

int SVD_svd_f (SVD me, float **m)
{
	long i, j;
	for (i = 1; i <= my numberOfRows; i++)
	{
		for (j = 1; j <= my numberOfColumns; j++) my u[i][j] = m[i][j];
	}
	return NUMsvdcmp (my u, my numberOfRows, my numberOfColumns, my d, my v);
}

void SVD_setTolerance (SVD me, double tolerance)
{
	my tolerance = tolerance;
}

double SVD_getTolerance (SVD me)
{
	return my tolerance;
}

int SVD_compute (SVD me)
{
	return NUMsvdcmp (my u, my numberOfRows, my numberOfColumns, my d, my v);
}

int SVD_solve (SVD me, double b[], double x[])
{
	if (! SVD_compute (me)) return 0;
	(void) SVD_zeroSmallSingularValues (me, 0);
	return NUMsvbksb (my u, my d, my v, my numberOfRows, my numberOfColumns, b, x);
}

int SVD_sort (SVD me)
{
	SVD thee = NULL; long i, j, *index = NULL;
	
	if (((thee = Data_copy (me)) == NULL) ||
		((index = NUMlvector (1, my numberOfColumns)) == NULL) ||
		! NUMindexx_d (my d, my numberOfColumns, index)) goto end;
			
	for (j = 1; j <= my numberOfColumns; j++)
	{
		long from = index[my numberOfColumns - j + 1];
		my d[j] = thy d[from];
		for (i = 1; i <= my numberOfRows; i++) my u[i][j] = thy u[i][from];
		for (i = 1; i <= my numberOfColumns; i++) my v[i][j] = thy v[i][from];
	}
end:
	forget (thee);
	NUMlvector_free (index, 1);
	return ! Melder_hasError ();
}

long SVD_zeroSmallSingularValues (SVD me, double tolerance)
{
	long i, numberOfZeroed = 0; 
	double dmax = my d[1];
	
	if (tolerance == 0) tolerance = my tolerance;
	for (i = 2; i <= my numberOfColumns; i++)
	{
		if (my d[i] > dmax) dmax = my d[i];
	}
	for (i = 1; i <= my numberOfColumns; i++)
	{
		if (my d[i] < dmax * tolerance)
		{
			my d[i] = 0; numberOfZeroed++;
		}
	}
	return numberOfZeroed;
}

/* End of file NUMclasses.c */
