litter

Doolitter分解三对⾓矩阵分解拟三对⾓分解

1#include

2#include

3#include

4#include

5#include

6#include

7#include

8

9

10classMclVector

11{

12public:

13intn;

14double*Mat;

15/**

16type=0:列向量(默认)

17type=1:⾏向量

18**/

19inttype;

20

21MclVector(){Mat=NULL;n=0;}

22MclVector(intlen,doubleinitVal=0.0)

23{

24n=len;

25Mat=newdouble[n+1];

26for(inti=0;i<=n;i++)Mat[i]=initVal;

27type=0;

28}

29

30

31doubleoperator[](intid)const

32{

33returnMat[id];

34}

35

36

37double&operator[](intid)

38{

39returnMat[id];

40}

41

42doublelength()const

43{

44doublesum=0;

45for(inti=1;i<=n;i++)sum+=Mat[i]*Mat[i];

46returnsqrt(sum);

47}

48

49MclVectoroperator*(doubleval)const

50{

51MclVectorans=MclVector(n);

52for(inti=1;i<=n;i++)ans[i]=Mat[i]*val;

53returnans;

54}

55

56MclVectoroperator/(doubleval)const

57{

58MclVectorans=MclVector(n);

59for(inti=1;i<=n;i++)ans[i]=Mat[i]/val;

60returnans;

61}

62

63MclVectoroperator+(constMclVector&newVector)const

63MclVectoroperator+(constMclVector&newVector)const

64{

65MclVectorans=MclVector(n);

66for(inti=1;i<=n;i++)ans[i]=Mat[i]+newVector[i];

67returnans;

68}

69

70

71MclVectoroperator-(constMclVector&newVector)const

72{

73MclVectorans=MclVector(n);

74for(inti=1;i<=n;i++)ans[i]=Mat[i]-newVector[i];

75returnans;

76}

77

78MclVectoroperator*=(doubleval)

79{

80for(inti=1;i<=n;i++)Mat[i]=Mat[i]*val;

81return*this;

82}

83

84MclVectoroperator/=(doubleval)

85{

86for(inti=1;i<=n;i++)Mat[i]=Mat[i]/val;

87return*this;

88}

89

90

91MclVectoroperator+=(constMclVector&newVector)

92{

93for(inti=1;i<=n;i++)Mat[i]+=newVector[i];

94return*this;

95}

96

97MclVectoroperator-=(constMclVector&newVector)

98{

99for(inti=1;i<=n;i++)Mat[i]-=newVector[i];

100return*this;

101}

102

103MclVectorGetTranspo()const

104{

105MclVectorans=*this;

=1;

107returnans;

108}

109

110

111voidprint()const

112{

113for(inti=1;i<=n;i++)printf("%8.3lf",Mat[i]);

114puts("");

115}

116

117

118};

119

120classMclMatrix

121{

122public:

123introw,col;

124MclVector*Mat;

125

126MclMatrix(){Mat=NULL;}

127MclMatrix(int_row,int_col,doubleinitVal=0.0)

128{

128{

129row=_row;

130col=_col;

131Mat=newMclVector[row+1];

132for(inti=0;i<=row;i++)Mat[i]=MclVector(col,initVal);

133}

134

135voidtIdentityMatrix()

136{

137for(inti=1;i<=row;i++)

138{

139for(intj=1;j<=col;j++)

140{

141if(i==j)Mat[i][j]=1;

142elMat[i][j]=0;

143}

144}

145}

146

147MclMatrixGetTranspo()const

148{

149MclMatrixans=MclMatrix(col,row);

150for(inti=1;i<=;i++)

151{

152for(intj=1;j<=;j++)

153{

154ans[i][j]=Mat[j][i];

155}

156}

157returnans;

158}

159

160voidprint()const

161{

162for(inti=1;i<=row;i++)Mat[i].print();

163puts("");

164}

165

166MclVector&operator[](intid)const

167{

168returnMat[id];

169}

170

171

172MclVector&operator[](intid)

173{

174returnMat[id];

175}

176

177MclMatrixoperator*(constMclMatrix&Matrix)const

178{

179MclMatrixans=MclMatrix(row,);

180for(inti=1;i<=row;i++)

181{

182for(intj=1;j<=;j++)

183{

184for(intk=1;k<=col;k++)

185{

186ans[i][j]+=Mat[i][k]*Matrix[k][j];

187}

188}

189}

190returnans;

191}

192

193MclMatrixoperator+(constMclMatrix&Matrix)const

193MclMatrixoperator+(constMclMatrix&Matrix)const

194{

195MclMatrixans=MclMatrix(row,);

196for(inti=1;i<=row;i++)

197{

198for(intj=1;j<=;j++)

199{

200ans[i][j]=Mat[i][j]+Matrix[i][j];

201}

202}

203returnans;

204}

205

206MclMatrixoperator-(constMclMatrix&Matrix)const

207{

208MclMatrixans=MclMatrix(row,);

209for(inti=1;i<=row;i++)

210{

211for(intj=1;j<=;j++)

212{

213ans[i][j]=Mat[i][j]-Matrix[i][j];

214}

215}

216returnans;

217}

218

219MclVectorGetCol(intcolId)const

220{

221MclVectorans=MclVector(row);

222for(inti=1;i<=row;i++)ans[i]=Mat[i][colId];

223returnans;

224}

225MclVectorGetRow(introwId)const

226{

227MclVectorans=MclVector(row);

228for(inti=1;i<=col;i++)ans[i]=Mat[rowId][i];

229returnans;

230}

231

232

233MclMatrixoperator*=(constMclMatrix&Matrix)

234{

235return*this=*this*Matrix;

236}

237MclMatrixoperator+=(constMclMatrix&Matrix)

238{

239return*this=*this+Matrix;

240}

241MclMatrixoperator-=(constMclMatrix&Matrix)

242{

243return*this=*this-Matrix;

244}

245

246MclMatrixoperator*(doublex)const

247{

248MclMatrixans=*this;

249for(inti=1;i<=row;i++)

250{

251for(intj=1;j<=col;j++)

252{

253ans[i][j]*=x;

254}

255}

256returnans;

257}

258

258

259};

260

261MclMatrixvectorMulVector(constMclVector&A,constMclVector&B)

262{

263if(==0)

264{

265MclMatrixans=MclMatrix(A.n,B.n);

266for(inti=1;i<=A.n;i++)

267{

268for(intj=1;j<=B.n;j++)

269{

270ans[i][j]+=A[i]*B[j];

271}

272}

273returnans;

274}

275el

276{

277asrt(A.n==B.n);

278MclMatrixans=MclMatrix(1,1);

279for(inti=1;i<=A.n;i++)

280{

281ans[1][1]+=A[i]*B[i];

282}

283returnans;

284}

285}

286

287intsgn(doublex)

288{

289if(x<-0.000001)return-1;

290if(x>0.000001)return1;

291return0;

292}

293

294

295

296/**

297矩阵的Doolittle分解:

298[1]Mat是⽅阵

299[2]Mat的前n-1阶主⼦式⾏列式不为0

300[3]分解的L为单位下三⾓阵

301[4]分解的U为上三⾓阵

302[5]返回值为

303**/

304std::pairDoolittleSplit(constMclMatrix&Mat)

305{

306intn=;

307MclMatrixL=MclMatrix(n,n);

308MclMatrixU=MclMatrix(n,n);

309for(intk=1;k<=n;k++)

310{

311for(intj=k;j<=n;j++)

312{

313U[k][j]=Mat[k][j];

314for(intt=1;t<=k-1;t++)U[k][j]-=L[k][t]*U[t][j];

315}

316if(k==n)continue;

317

318for(inti=k+1;i<=n;i++)

319{

320L[i][k]=Mat[i][k];

321for(intt=1;t<=k-1;t++)L[i][k]-=L[i][t]*U[t][k];

322L[i][k]/=U[k][k];

323}

323}

324}

325for(inti=1;i<=n;i++)L[i][i]=1;

326returnstd::make_pair(L,U);

327}

328

329

330

331/**

332三⾓矩阵分解:

333[1]Mat是⽅阵

334[2]jr时Mat[i][j]=0

335[2]j>i且j-i>s时Mat[i][j]=0

336**/

337std::pairTriangleSplit(constMclMatrix&Mat,intr,ints)

338{

339intn=;

340MclMatrixL=MclMatrix(n,n);

341MclMatrixU=MclMatrix(n,n);

342for(intk=1;k<=n;k++)

343{

344for(intj=k;j<=n;j++)

345{

346U[k][j]=Mat[k][j];

347for(intt=std::max(1,std::max(k-r,j-s));t<=k-1;t++)U[k][j]-=L[k][t]*U[t][j];

348}

349if(k==n)continue;

350

351for(inti=k+1;i<=n;i++)

352{

353L[i][k]=Mat[i][k];

354for(intt=std::max(1,std::max(i-r,k-s));t<=k-1;t++)L[i][k]-=L[i][t]*U[t][k];

355L[i][k]/=U[k][k];

356}

357}

358for(inti=1;i<=n;i++)L[i][i]=1;

359returnstd::make_pair(L,U);

360}

361

362/**

363拟三对⾓矩阵分解

364对n=5矩阵A样⼦如下：

365a1c100d1

366d2a2c200

3670d3a3c30

36800d4a4c4

369c500d5a5

370即输⼊为三个长度为n的向量

371

372A=LU

373L样⼦如下：

374p10000

375d2p2000

3760d3p300

37700d3p40

378r1r2r3r4r5

379

380U样⼦如下:

3811q100s1

38201q20s2

383001q3s3

3840001s4

38500001

386

387即将返回p,q,s,r四个向量(所有的向量长度都是n)

388vector[0]=p

388vector[0]=p

389vector[1]=q

390vector[2]=s

391vector[3]=r

392**/

393

394std::vectorQuasiDiagonalSplit(constMclVector&a,constMclVector&c,constMclVector&d)

395{

396intn=a.n;

397asrt(c.n==n);

398asrt(d.n==n);

399asrt(n>0);

400

401MclVectorp=MclVector(n);

402MclVectorq=MclVector(n);

403MclVectors=MclVector(n);

404MclVectorr=MclVector(n);

405

406p[1]=a[1];

407for(inti=1;i<=n-2;i++)

408{

409q[i]=c[i]/p[i];

410p[i+1]=a[i+1]-d[i+1]*q[i];

411}

412

413s[1]=d[1]/p[1];

414for(inti=2;i<=n-2;i++)s[i]=-d[i]*s[i-1]/p[i];

415s[n-1]=(c[n-1]-d[n-1]*s[n-2])/p[n-1];

416

417r[1]=c[n];

418for(intj=2;j<=n-2;j++)r[j]=-r[j-1]*q[j-1];

419r[n-1]=d[n]-r[n-2]*q[n-2];

420r[n]=a[n];

421for(intj=1;j<=n-1;j++)r[n]=r[n]-r[j]*s[j];

422

423std::vectorans;

_back(p);

_back(q);

_back(s);

_back(r);

428

429returnans;

430}