SVD (c语言代码)

/*
* svdcomp - SVD decomposition routine.
* Takes an mxn matrix a and decomposes it into udv, where u,v are
* left and right orthogonal transformation matrices, and d is a
* diagonal matrix of singular values.
*
* This routine is adapted from svdecomp.c in XLISP-STAT 2.1 which is
* code from Numerical Recipes adapted by Luke Tierney and David Betz.
*
* Input to dsvd is as follows:
*   a = mxn matrix to be decomposed, gets overwritten with u
*   m = row dimension of a
*   n = column dimension of a
*   w = returns the vector of singular values of a
*   v = returns the right orthogonal transformation matrix
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define SIGN(a,b) ((b) >= 0.0 ? fabs(a) : -fabs(a))
#define MAX(a,b) (a>b?a:b)
#define MIN(a,b) (a<b?a:b)
float **Make2DArray(int row, int col){
	float **a;
	int i;
	a = (float **)malloc(row * sizeof(float *));
	for (i = 0; i<row; i++){
		a[i] = (float *)malloc(col * sizeof(float));
	}
	return a;
}
void SetZero(float** a, int m){
	int i, j;
	for (i = 0; i < m; i++){
		for (j = 0; j < m; j++){
			a[i][j] = 0;
		}
	}
}
void Free2DArray(float** a, int row){
	int i;
	for (i = 0; i<row; i++){
		free(a[i]);
	}
	free(a);
}

void p(float**a, int m, int n){
	for (int i = 0; i < m; i++){
		for (int j = 0; j < n; j++)
			printf("%f ", a[i][j]);
		printf("\n");
	}
}

static double PYTHAG(double a, double b)
{
	double at = fabs(a), bt = fabs(b), ct, result;

	if (at > bt)       { ct = bt / at; result = at * sqrt(1.0 + ct * ct); }
	else if (bt > 0.0) { ct = at / bt; result = bt * sqrt(1.0 + ct * ct); }
	else result = 0.0;
	return(result);
}

int dsvd(float **a, int m, int n, float *w, float **v)
{
	int flag, i, its, j, jj, k, l, nm;
	double c, f, h, s, x, y, z;
	double anorm = 0.0, g = 0.0, scale = 0.0;
	double *rv1;

	if (m < n)
	{
		fprintf(stderr, "#rows must be &g