#include <math.h>
#include "mex.h"
double mymax(double x, double y)
{
if (x > y)
return x;
else
return y;
}
double absolute(double x)
{
if (x >= -x)
return x;
else
return -x;
}
void permuteInt(int *x, int p, int q)
{
int temp;
temp = x[p];
x[p] = x[q];
x[q] = temp;
}
void permute(double *x, int p, int q)
{
double temp;
temp = x[p];
x[p] = x[q];
x[q] = temp;
}
void permuteRows(double *x, int p, int q,int n)
{
int i;
double temp;
for(i = 0; i < n; i++)
{
temp = x[p+i*n];
x[p+i*n] = x[q+i*n];
x[q+i*n] = temp;
}
}
void permuteCols(double *x, int p, int q,int n)
{
int i;
double temp;
for(i = 0; i < n; i++)
{
temp = x[i+p*n];
x[i+p*n] = x[i+q*n];
x[i+q*n] = temp;
}
}
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
int n,sizL[2],sizD[2],i,j,q,s,
*P;
double mu,gamma,xi,delta,beta,maxVal,theta,
*c, *H, *L, *D, *A;
/* Input */
H = mxGetPr(prhs[0]);
if (nrhs == 1)
{
mu = 1e-12;
}
else
{
mu = mxGetScalar(prhs[1]);
}
/* Compute Sizes */
n = mxGetDimensions(prhs[0])[0];
/* Form Output */
sizL[0] = n;
sizL[1] = n;
plhs[0] = mxCreateNumericArray(2,sizL,mxDOUBLE_CLASS,mxREAL);
L = mxGetPr(plhs[0]);
sizD[0] = n;
sizD[1] = 1;
plhs[1] = mxCreateNumericArray(2,sizD,mxDOUBLE_CLASS,mxREAL);
D = mxGetPr(plhs[1]);
plhs[2] = mxCreateNumericArray(2,sizD,mxINT32_CLASS,mxREAL);
P = (int*)mxGetData(plhs[2]);
/* Initialize */
c = mxCalloc(n*n,sizeof(double));
A = mxCalloc(n*n,sizeof(double));
for (i = 0; i < n; i++)
{
P[i] = i;
for (j = 0;j < n; j++)
{
A[i+n*j] = H[i+n*j];
}
}
gamma = 0;
for (i = 0; i < n; i++)
{
L[i+n*i] = 1;
c[i+n*i] = A[i+n*i];
}
/* Compute modification parameters */
gamma = -1;
xi = -1;
for (i = 0; i < n; i++)
{
gamma = mymax(gamma,absolute(A[i+n*i]));
for (j = 0;j < n; j++)
{
//printf("A(%d,%d) = %f, %f\n",i,j,A[i+n*j],absolute(A[i+n*j]));
if (i != j)
xi = mymax(xi,absolute(A[i+n*j]));
}
}
delta = mu*mymax(gamma+xi,1);
if (n > 1)
{
beta = sqrt(mymax(gamma,mymax(mu,xi/sqrt(n*n-1))));
}
else
{
beta = sqrt(mymax(gamma,mu));
}
for (j = 0; j < n; j++)
{
/* Find q that results in Best Permutation with j */
maxVal = -1;
q = 0;
for(i = j; i < n; i++)
{
if (absolute(c[i+n*i]) > maxVal)
{
maxVal = mymax(maxVal,absolute(c[i+n*i]));
q = i;
}
}
/* Permute D,c,L,A,P */
permute(D,j,q);
permuteInt(P,j,q);
permuteRows(c,j,q,n);
permuteCols(c,j,q,n);
permuteRows(L,j,q,n);
permuteCols(L,j,q,n);
permuteRows(A,j,q,n);
permuteCols(A,j,q,n);
for(s = 0; s <= j-1; s++)
L[j+n*s] = c[j+n*s]/D[s];
for(i = j+1; i < n; i++)
{
c[i+j*n] = A[i+j*n];
for(s = 0; s <= j-1; s++)
{
c[i+j*n] -= L[j+n*s]*c[i+n*s];
}
}
theta = 0;
if (j < n-1)
{
for(i = j+1;i < n; i++)
theta = mymax(theta,absolute(c[i+n*j]));
}
D[j] = mymax(absolute(c[j+n*j]),mymax(delta,theta*theta/(beta*beta)));
if (j < n-1)
{
for(i = j+1; i < n; i++)
{
c[i+n*i] = c[i+n*i] - c[i+n*j]*c[i+n*j]/D[j];
}
}
}
for(i = 0; i < n; i++)
P[i]++;
mxFree(c);
mxFree(A);
}
IT狂飙
- 粉丝: 4842
- 资源: 2650
最新资源
- (2025)国家基层糖尿病防治管理指南认证考试试题及答案.docx
- (2025)国家公务员录用考试行测常识题库及答案.docx
- (2025)汉字听写大会试题库(附答案).docx
- (2025)国家开放大学《中国法律史》形成性考核1-4与参考答案.docx
- (2025)工业机器人技术题库及答案.docx
- (2025)科创板股票投资知识题库及答案.docx
- (2025)护理三基基础知识考试题库(含答案).docx
- 知识领域:仪器仪表,变流器,自动控制 关键词:光伏MPPT,电压控制器,微电网,河南求同电气,光伏模拟实验系统
- STM32F107各种接口程序合集工程文件 包含串口,CAN,时钟芯片,FLASH,外包AT25320储存,数据结构,枚举,适合刚刚出来工作的工程师以及进阶工程师 1.提供AD STM32F107原
- 基于fpga的多功能pwm模块设计 可应用于:dab,llc,buck,boost,全桥,推娩等dcdc电路 功能: 1.输出多路互补的pwm 2.每路互补pwm死区可调 3.每路互补pwm频率独立
- nianhuishougao
- 三菱PLC分拣程序基于三菱FX系列的分拣程序,可用于学习
- 电力系统的物理信息神经网络python源代码 代码按照高水平文章复现 介绍了一种在电力系统中应用物理信息神经网络的框架 利用控制电力系统的基本物理定律,并受到机器学习领域最新发展的启发,我们提出了一
- 知识领域:变流器,自动控制 关键词:软锁相环,河南求同电气,电压不平衡,微电网并网系统,变流器
- 三菱FX3U与4台英威腾GD系列变频器通讯案例实战程序 有注释,并附送程序,有接线方式,设置 器件:三菱FX3U的PLC,4台英威腾GD系列变频器,昆仑通态 功能:实现频率设
- 中颖正弦波矢量电动车控制器 1-提供原理图 2-提供pcb图 3-提供C源代码(主芯片SH79F3213) 带自学习功能,可任意匹配电机
资源上传下载、课程学习等过程中有任何疑问或建议,欢迎提出宝贵意见哦~我们会及时处理!
点击此处反馈