哈尔滨某大学结构力学I一一

SubroutineElem_Stiff(······)

说明Stiff=0.0!单元刚度清零SelectCase(Type)Case(1)

平面杆系结构单元Case(2)

空间杆系结构单元CaseDefault

出错信息EndSelectEndSubroutineElem_Stiff3.4杆系结构单元分析子程序3.4.1单元刚度总体设计3.4.2说明部分设计Integer,Intent(in)::···

入口整型参数Real(8),Intent(in)::···

入口实型参数Real(8),Intent(out)::···

出口实型参数Real(8)::Work1,···Integer::i,j,k,···实型和整型工作变量3.4杆系结构单元分析子程序3.4.3平面杆系结构设计SelectCase(Plane)Case(1)

平面桁架元素赋值Case(2)

平面梁柱元素赋值Case(3)······CaseDefault

出错信息EndSelect3.4.4空间杆系结构设计SelectCase(Space)Case(1)

空间桁架元素赋值Case(2)

空间梁柱元素赋值Case(3)

交叉梁元素赋值CaseDefault

出错信息EndSelect3.4杆系结构单元分析子程序3.4.5有关单元等效结点荷载设计和进一步的考虑1)单元等效结点荷载设计同仿单元刚度。2)从各类单元刚度元素的计算，可看到要用到长度、单元弹性特性、单元截面特性等数据。因此，要确定存放它们的数据结构。要将它们作为出口。3)为计算单元等效结点荷载元素，首先要建立各种荷载情况等效荷载表达式，它们可由积分或载常数表得到。然后要解决荷载信息的存放结构，也要将它们作为出口量。4)单元刚度矩阵、等效结点荷载矩阵都应先清零。4.1杆系结构整体分析

首先就全刚结点平面刚架进行讨论，然后推广。4.1.1总的思路在单元特性搞清后，将单元拼装回去。在结点处位移自动协调基础上，如果全部结点平衡，则求得的结点位移将是实际结构的解。因此，整体分析就是设法建立结点平衡方程。4.1.2坐标转换组成结构的杆件可以各个方向，单元分析对局部坐标，因此，必须将物理量转为统一坐标——整体坐标。1)力的转换关系4.1杆系结构整体分析2)位移转换关系3)转换矩阵转换矩阵是正交矩阵。4.1杆系结构整体分析4)杆端力转换5)杆端位移转换6)刚度方程的转换如果记称为整体单元刚度矩阵则这就是整体坐标下的单元刚度方程。本节以后的讨论认为都是对整体坐标的4.1杆系结构整体分析4.1.3结点平衡方程的建立1)一简单例子（如图）图中有两套编号，红的是单元杆端编号，黑的是结构整体编号。1-1)结点示意121221①②③图中蓝色的表示结点荷载（已知），红色的表示杆端力（未知的），、分别1、2单元杆端力子矩阵。对1、4结点“荷载”含有未知反力。21-2)结点平衡由示意图可见，结构结点的平衡方程为4.1杆系结构整体分析从例图可见，其全部结点平衡方程为121221①②③若记24.1杆系结构整体分析式中[I]、[0]分别为单位和零矩阵。若引入矩阵记号，则结点平衡方程可改写作这一结论虽然是由一个例子得到的，但是显然对一切结构都是成立的。问题在于不同结构，[A]矩阵是不同的。4.1杆系结构整体分析4.1.4杆端位移用结点位移来表示121221①②③仍以简单例子来说明若记由结点、杆端位移的协调条件，可得[]、[]的对应关系为式中[A]T是前面力关系[A]的转置，因此[A]T称为位移转换矩阵。4.1杆系结构整体分析4.1.5整体刚度方程——结点平衡121221①②③若记引入位移转换关系，则这就是整体刚度方程，它的物理实质是结点平衡。[K]称作结构刚度矩阵（或整体刚度矩阵），[R]称作综合等效结点荷载矩阵，它由两部分组成。单元个数4.1杆系结构整体分析4.1.6整体刚度矩阵的建立121221①②③若将[A]按单元分成图示三个子矩阵则由此可见，整体刚度矩阵可由各单元整体刚度矩阵装配累加得到。为说明如何装配，先将单元刚度矩阵进行分割整体结点码则由矩阵乘法可证明，[A]i[k]i[A]iT的结果是，将刚度矩阵子矩阵按整体结点码r、s送整体刚度矩阵相应位置。这一装配规则称为“对号入座”。4.1杆系结构整体分析1)任意结构情况上面结论是通过具体例子（全刚结点平面刚架）得到的，由虚位移原理或势能原理进行整体分析（见讲义），可得任意结构其结论同此例。2)结点位移编号如果按结点顺序，对结点非零位移进行依次编号，这一序号称作结点位移码。为便于计算机处理并减少结构刚度矩阵的阶次，将零位移的号码变为零。对图示三铰刚架，当仅用一种单元（梁柱自由是单元）时结点位移编号如图所示。3)单元定位向量按单元局部结点码顺序，将结点位移码排成的向量，称作单元的定位向量。①②③④①②③④4.1杆系结构整体分析对图示刚架各单元的定位向量为①(0,0,1,3,4,5)②(0,0,2,10,11,12)③(3,4,5,6,7,8)④(6,7,9,10,11,12)①(0,0,1,2,3)②(0,0,,6,7,8)③(1,2,3,4,5)④(4,5,6,7,8)4)按单元定位向量集装刚度矩阵和综合荷载前面说明的是分块子矩阵集装,下面说明如何按定为向量来集装.如果如图所是采用各种不同的单元（一端有铰），则定位向量为①②③④①②③④如何获得带铰的单元刚度矩阵和等效荷载矩阵定位向量①②③④①②③④4.1杆系结构整体分析4-1)刚度集装以3单元为例来说明定位向量单元局部位移码根据单元局部位移码和定位向量的对应关系用定位向量位移码送元素。根据单元局部位移码和定位向量的对应关系用定位向量位移码送元素，定位向量元素为零时不送。①②③④①②③④4.1杆系结构整体分析4-2)荷载集装以4单元为例来说明定位向量局部位移码此结论同样适用于刚度集装4.1杆系结构整体分析4.1.7整体分析总结1)对局部坐标和整体坐标不一致的单元，要对刚度、荷载进行坐标转换。2)需对“结构”进行结点、位移的局部和整体编号。4)集装所得整体刚度矩阵是对称、带状稀疏矩阵，当支撑条件能限制刚体位移时，矩阵非奇异。3)根据单元局部位移码和定位向量的对应关系用定位向量位移码送元素，定位向量元素为零时不送。据此可集装、累加得到整体刚度矩阵。5)综合荷载由两部分组成，因此首先要将直接作用结点的荷载按结点位移码送入，如果还有单元等效荷载，再按定位向量集装、累加。4.1杆系结构整体分析8)如果有某位移码方向弹性支撑，需进行将弹簧刚度送入位移码对应的对角线元素位置累加。9)如果有某位移码方向已知支撑位移，需进行将“边界条件处理”。具体做法以后介绍。7)整体刚度方程实质是全部结点的平衡条件。6)刚度矩阵带状稀疏，其带宽取决于结点、位移编码。最大半带宽=定位向量中最大元素差+1。4.1杆系结构整体分析2.5.8边界条件的处理10)当用虚位移或势能原理作整体分析时（势能为例），应变能为单元应变能之和，外力势能为单元外力势能之和+结点外力势能。全部杆端力的外力势能彼此抵消。1)乘大数法设第i个位移为已知值a，N=108或更大的数。乘大数法是将刚度矩阵Kii改为NKii，将Ri改为

Na。请考虑为什麽这样做能使边界条件得到满足？2)置换法（划零置1）设第i个位移为已知值a。4.1杆系结构整体分析上述置换工作量大一些，显然可看出边界条件得到精确满足。4.1杆系结社构整体宅分析★3)关于斜丛边界的短处理如图贤示意都的斜扛支座耗情况主，有却多种御处理孝方案亡。3-1撑)通过单具元的坐材标转换睛来处理xy图示赔有斜偿支座懂单元劣，r结点窝处以倾孔角-来进跌行坐顷标转元换，纽奉也即在r结点处皆整体坐制标为图斑示xy。r3-遗2)通过增鸡加一个洽单元来雾处理图示有槽斜支座伴单元，r结点处补沿y方向把增加强一个刚结的枕单元，此单元有“无炎穷大绝”的江抗拉红刚度撒、但没有抗类弯刚度适。单元战长度可艰任意。3-鸦3)对整体浮刚度矩丑阵进行薄处理（认参见匡采文起教薪材）最大半带宽2.5杆系结轨构整体吗分析2.抽5.巩10总刚度听矩阵的悦存储与踪蝶对应解回法因为邻总刚傍度矩惜阵对拍称、吨带状疾稀疏梳，利附用这岁一特点可减今少存储苗、提高孤解算效队率。零元素零元脱素非零元素最大半带宽主对角线元素等半带存储零元素非零元素变带宽一维存储带宽是变的到P牌.5祥52.5杆系支结构读整体业分析目前一腰般都用单变带宽砍存储，蜡下面结添合程序筹说明存储引和解偿法。咐首先酱介绍呼一些F9干0的语法稻。定义进导出型类型导出类辰型——结点typ光e:缎:typ洪_Jo辉intrea旁l色::口x,y！坐标in槽te坝ge铅r亦::尸G服DO荡F(死3)！整体混位移忌码end骄ty恼pety传p_洋Jo殊in衫t1)有关F90语法导出类摄型新特愤性用结禁点导哭出类肥型作州为成奏员导勒出单近元类爪型:typ厘e:忙:typ社_El表eme算ntint婶ege观r陶::输Joi膝ntN副o(2够)！结点编姥号ty颈pe郊(t眨yp冰_J搜oi极nt)牵::右No孟de(案2)！结点int蚀ege逝r范::滴Glb午DOF黎(6)认!定位盟向量re污al弄::秩E贿A,绞EIend怎ty仙pety始p_丙El达em亩en妙ttyp掏e(ty晨p_铺El助em趋en汉t)过::泻Ele亡m(5及)！定义5个单斥元类休型…!对单良元i的端孤点j的x,y凳,GD区OF(译1:3拦)的赋值Ele份m(i舰)%N庭ode潜(j)=Joi以nt(腔Ele踪蝶m(i缝)%J膝oin废tNo碰(j))…由导出沉类型定百义新类葬型由导出白类型定肃义变量re肺al冲:蜡:互A(笨5)咏,B福(5口,1辟0)水,C股(5从)B=0窝.0!对B清零A=虾1.汽0!对A赋1:勾A出(i琴)=般1.料0,伙i绞=1引,5C=圈A+招2!数组貌与标惯量运误算:A堤(1:执5)+忠(/抢2,2合,2,企2,2柴/)A=纲C+贪A!数组昆与数亚组运等算（鹿同形贴）C=sq证rt喝(A)!数组径的函捉数运殃算:C寇(i)赠=sq醒rt颂(A稠(i),颂i=1冲,5数组内尘部函数：do矛t_谊pr摸od程uc床t(存ve谣ct全or圾_a仅,v萄ec胆to夫r_卖b)!点积如：dot召_pr乡丰odu航ct(凤(/1衬,2,穗3/)沾,(/焦2,3甚,4/妹))的值为20（待价续）有关F9水0语法数组运真算与赋业值：ma塞tm恩ul纠(m红at坏ri弓x_姑a,吹ma碰tr刘ix睁_b)矩阵次相乘如：lo集cE盟Di滨sp=mat侨mul僚(T,度glb粱EDi矩sp)tra兔nsp落ose船(ma告tri拣x)矩阵转泉置如：glb犁EDi浅sp=ma食tm储ul虎(t择ra冠ns占po鼓se浙(T亮),晚lo好cE萌Di极sp)siz矿e(a死rra个y,d洪im)求数组桥第di俊m维的汉长度dim为可昌选变框元:si内ze捷(a住,d糟im兵=2调)若arr绘ay为一维栽时，可小不用di姻m。sum糠(ar叛ray袜,di初m,m席ask嘴)数组元扬素求和di苏m,匠ma岔sk为可盐选变侵元；mas欧k=条件表犬达式sum软(a(锣1:1爽0))对a的1到10元素求容和su军m(午a(突1:淡10悼),阶ma伶sk追=a折>0摸)对a(1剂:10台)中大于0的元咬素求包和（续）有关F9见0语法whe位re结构新特铁性例whe幅re乓(C>科0)C=灿0A=柴B*贼Dend目wh臭erewhe摸re轰(C>焰0)A恢=坛Ben检d迹wh辣er搞e定义wh渗er禾e(数组关六系表达防式)数组齿赋值姨语句…els逆ewh镇ere数组泛赋值社语句...end苍wh驾ere规则:1)同形数犁组;2摆)不许环嵌套;3材)最多二洗分叉有关F9嫂0语法cy肯cl矛e和exi汗t语句新特或性用在do循环穗中cy丧cl肢e——作下一规个循环收步exi庙t——跳出涝循环坟，执有行en锣d接do后一哈条语眨句等效例do...if倒(.洲no停t.烂co艳nd寄1)区th覆en...if(跨con本d2)got宏o5..止.end舟ifen递d售do5柜..巩.用法do..遥.if骨(咸co仙nd孕1)票c头yc透le...if驰(co克nd2弯)e甩xit...end探do...有关F9大0语法数组捏构造汉函数spr您ead语法spr杰ead歪(数组名,di巡寿m,吼nc萄op叨ie该s)将数组删沿dim维方丘向复惭制nc吉op咽ie异s形成新张数组di律m,粗nc艇op死ie摔s—整型、暑位置变倚元、关掀键字变斗元（若无按位衫置引津用，汽可略奸关键秩字）例:（仅辉限一其维数盒组）1）spr输ead围(on焰e,d君im=泄1,n窝cop舰ies铁=3)spr锐ead葬(on草e,1粒,3)spr救ead浴(on颗e,n羊cop牺ies=3,赵dim嚷=1)[1,1,1]或[1,1,1]T2）ELo讨cVe盘c(1跨:6)舟=(/浆1,脑0,3绞,4,规5,0裂/)spr暮ead懂(EL左ocV低ec,逐dim=1,串nco察pie皇s=3举)3）spr惜ead众(A(吊2,2拳:),取dim具=1,垂nco厘pie梨s=2叔)如果di燥m=景2呢?有关F9丸0语法指针po慨in糖te冈rpo趁in叉te搭r是变量搬的属性眼，可以旬指向相社同类型谜的变量;被指向叙的目标块必须具寄有ta袍rg明et属性或poi乎nte课r属性可以将话指针变注量理解半为别名赔、称号re州al闲,t锯ar经ge考t棒:李:肚a,并b,殃ED宪is焰p(巩6)!可被略指针蜘所指re竿al视,p厕oi怪nt陶er诱:秃:养p1厌,p侄2!称号补：班彼长、芝课代联表！p1丘,p物2是指针屿，可以采指向实课型数据re永al吓,p含oi造nt夹er堪:朴:他G(晋:)!先进集承体！G是指君针，少可以舍指向落一维开实型论数组指针融是一酿种“称号”，上丈述声抗明语珠句建制立了“称号骄”，但并培未“授意予”某个呜变量续这个性称号貌，因但此是典指向“空”，并未增占用内仅存a=性3.赖0p1牙=>哭a!p身1指向a!称号p1授予a，a的数据盾有两个功名：固示定名a和流动缘瑞名p1逃;既可用p1也可盘用a（p1丙—班长，a立—张三）a=饰4.榨0!者a的值变馆为4.咬0,p1也变米为4.0p1编=职>朱b!班长换再人了G丑=>ED扭is鹊p!先进集渐体有了盐得主！ED芝is齿p(:)的长度停就是G(:饭)的长淋度，赚用G和用EDi睛sp同样效匹果又如：rea晃l,t若arg容et伐::慰a,brea菌l,p糕oin毁ter名::钳p,qa=指3.卫14b=骄2.愚0p究=>去a!丹p哨=沿a顿=讨3.工14q=瞎>b!沿q岔=秧b膛=鉴2.啊0q坝=横p!昼(q指向顶的数奔据b)=崇(p指向的每数据a梁)！即：岭所有=3箩.14指针集可以晃指向夺有名映的数态据区厌，也捆可以拼指向原无名堡的数戚据区rea舅l,p护oin肺ter筒::洗b(:迷,:)in杠te帖ge计r罩::新nre骗ad鼠(*粒,*拌)牲nal拜lo都ca画te法(津b(钞n,演n)启)!指针蔽指向闲了一逼个刚角开辟乳的数叙组!以下可泻以当作阳常规数福组用b(奸1,丧1)织=情1腿.1b(1机,2)泉=索1.2...dea廉llo脾cat拿e(b呜)有关F90语法用指封针建坚立动普态数倡组指针三与al摘lo像ca坚ta纳bl白e数组的描区别具备al管lo雄ca祸ta停bl丢e数组的咳所有功骄能还可以洪用在导到出类型株中，例汇如整体坦刚度矩弯阵的变挖带宽存境储：ty茫pe消:亲:ty绸p_泼Kc卫ol!整体刚角度矩阵K的列rea材l(8踏),p尿oin旦ter予::笋ro妙w(:烫)!该列俘的行判元素en夹d睬ty积pe..扩.typ蜡e(ty播p_若Kc刚ol间),惑al睬lo跪ca堵ta尼bl谣e::Kc脸ol(:)…all抓oca讽te纷(Kc扰ol真(N疼Gl巴bD刮OF))!分配了NG就lb暴DO活F列...all倍oca阶te滩(Kc浴ol(厌5)%萄row搏(3:彩5))!第5列只用3至5行(1奥)NE昼le吵m,裙NJ疾oi害nt,NG松lb计DO面F,NJL患oad,NE永Lo浪ad单元数画结点逼数总雨自由度贵数结凝点荷载肚数单箭元荷载雷数[Jo什int然-结点]…NJ雨oi误nt行(2)纲Jo尚int壤%X,摘Joi谜nt%那Y,这GDO介F(1盆:3)X坐标Y坐标睁结岭点位移旺码[El烘em蜡-单元]…NE肝le果m行(3)歌Jo醋int加No1设,Jo体int峰No2卷,E晨A,E偿I起点沃号外终点品号次刚度[JL俱oa鹅d-结点荷纸载]…NJ惕Lo植ad行(4)Jo交in易tN还o,Lod右DOF,Lod的Val作用点援号局怀部位移阴码荷翠载值[EL本oa胞d-单元荷哥载]…NE悟Lo竞ad行(5)El游em慢No,Ind拦x,印Po爸s,Lod忆Val单元号睡类型躁号位仅置荷勇载值Indx

类型

pos1均布长度2集中位置3...2)某程序莲输入数态据说明3,5爹,7,游1,10,0鞋,0即,0,董00,载4,算1深,2顷,34,4康,4践,5,的64,4半,4线,5,撞74,0挪,0夫,0,希01,2迁,1恳.0E迟9,季162,3萌,1曲.0E搁9,混245,拦4,友1哀.0蜂E9刚,防122,傻1,向10汇.0槐E3!结点2，自欧由度1，值为10E跌32,或1,扎4,哥-4彻.0察E3！单元2，均汇布，吊长4米，值大为-4馅E32-1尖)数据亲文件捐例子宰：(2)(1)(3)24135i=6i=4i=310kN4kN/m4m4mEA=恢109N(1编)(2)坐标位移宣码(3赠)(4)(5)结点烟号EA,曾EIrea薄d(5酿,*)NEl罚em,延NJo捐int尖,NG宴lbD墙OF,斧NJL饱oad砖,NE繁Loa肿dall及oca涛te建(Jo耗in叠t(径NJ回oi换nt))all庭oca朽te排(Ele触m(N巩Ele糖m))all数oca林te拳(JLo泰ad(域NJL询oad))all配oca贞te芹(EL肾oa帅d(经NE盆Lo萄ad))..录.rea惑d(5扯,*)蹄(J字oin匀t(i揉),i狠=1,舒NJo摸int歉)re剥ad脖(5绕,*畜)验(Ele欺m(i关e)%袜Joi针ntN撞o,E装lem厉(ie求)%E左A,&ELe垄m(i嗽e)%箩EI,出ie=1,晨NEl递em)if桐(NJ后Lo咳ad>0)羽re汗ad(狭5,*京)(JL伤oa性d(活i)棚,i=1宏,N矛JL堆oa农d)if巡寿(NE劲Lo捕ad>0)冒re岔ad(若5,*果)(ELo雁ad(为i),俘i=1,供NEL炒oad恒)2-酸2)程序读镰入计算病所需数跌据：2-3楼)求始行床码和分扯配带宽邪子程序!=春==紫==就==改==犹==坊==血==皮==随==慕==杰==援==故==壮==返==拾==贱=sub塑rou仁tin老eSet宫Mat津Ban寇d(Kc幸ol,E挣lem左)!接口里简单!=划==幼==鸽==鸭==岁==坐==进==前==烧==斑==姥==民==遵==断==伤==循==扛=typ致e(t牧yp_雕Kco麦l),延int哥ent曾(inou序t)缸:室:Kco判l(:)！总刚列ty湖pe扩(t劫yp款_E疏le躬me皆nt予),粱in持te裕nt掠(i晚n))湖:圆:滥El敏em叮(:霉)！单元int重ege沾r(i灯kin围d)夕::扒m检in羊DO样F,盏EL楼oc举Ve掠c(开6)in该te赖ge吨r(位ik于in盟d)询::诸R似ow耍1(菜si截ze哪(K输co缩慧l,继di远m=诸1)扣)！Row惧1为自练动数畜组，揉子程岭序结碗束后摘自动乒释放债。！这脉样做壤可使敬接口候简单快，减湾少了溜数组卖的控慈制变笼量。int助ege肯r(i闷kin妻d)溜::ie每,j饰,N阶Gl浩bD售OF做,N殃El哈emNG宿lb己DO抬F=siz滥e(K居col区,di初m=1)！使接口王简单NE扔le申m=siz日e(E序lem齿,di住m=1岔)Ro热w1踢=NG吩lb营DO植F!先设所科有始行滴码同终悲行码!确定各写列始行侨码向量doie=1列,N浆el盾em!对单元穗循环!确定摘定位陡向量EL陆oc缺Ve柔c(:龙)很=El此em守(i摩e)捷%G锡lb冰DO习F(:)!寻找虾定位网向量左中大品于零罢的最桂小值min婶DOF=min宿val手(EL索ocV钉ec,皮mas城k=ELo省cVe衰c>0膨)!屏蔽徒定位用向量乌中小庸于等薯于零泥的wh姻er被e隐(ELo点cVe父c>腔0)!寻找Ro番w1匆(E饿Lo吗cV拔ec多)和mi逆nD感OF中的最红小值Row盾1(E咬Loc等Vec粱)=苍mi宣n(R南ow1慕(EL谁ocV偏ec)苏,&mi揉nD钳OF)end槐wh浇ereen默d梁do!为各列创的带宽犹分配空武间do脊j看=1吓,N洽Gl槽bD比OF!对总自气由度数呀循环all味oca锅te酱(Kc淹ol(卵j)%忧row河(Ro薪w1(读j):腥j)!给Kc担ol佣(j事)%奶ro蚂w从Ro骄w1难(j绩)到j个空陕间Kco删l(j涉)%r巷ow=z市ero!总刚元田素清零end萝doret灶urnend勤su抵bro急uti炎neSet披Mat棒Ban娱d3)整体刚游度矩阵边的集成doie=1,不NEl椅em…!计算单阻刚EK苍(6劲,6种),沫EL钞oc蜘Ve刘c(产6)do反j=1看,6!对单元截逐列集岩成JGD耍OF凉=EL裕oc测Ve艘c(谈j)!取出位监移码if戒(JG言DOF神=丹0)飞cyc游le!作下一码循环步whe快re鄙(ELo厉cVe醋c>0卧.娃And翠.ELo远cVe亏c<=秃J悼GD增OF斩)!位移煮码非侨零同旁时小荒于第j个局孟部码欠对应爪的位顶移码Kco尚l(J雨GDO闹F)%殃row闸(EL迎ocV申ec)葬=构&Kc师ol墨(J剖GD轧OF隶)%灶ro驴w(蔬EL衡oc趋Ve卫c)答+你EK胖(:把,j幕)!集成en排d暂wh能er牺eend灰doen挖d隆do局部码位移码4)变带宽猴矩阵的张分解求锄解4-励1)LD田LT分解法艺求解[A]{x}={b}若对称正定，则可分解为=单位歌上三姨角阵对角傍阵原方惯程变肌为求解步政骤：1.分解：2.解y：3.解x：LD脸LT分解法Gau烤ss消去前消去处理缩慧右端昏项回代(不同的b只做弱一次蚂分解)存放粗：主对角—上三壁角—时不必后求和(上三怎角:i<j)不动存在处存在处(第j列系数)上三堂角:i<j4-育2)分解胶一般惯情况历：逐列忠分解对角线:i=j新>1(第j列系鼻数)第i列中第1个非岸零元坝素的嚼行码祸为：4-乖3)变列宽夜存贮的降分解修肾正：第j列中鲁第1个非零末元素的谦行码为封：则分解调顺序4-药4)仍F搁90实现!三角躲分解dia丝式g(1霞:nc甘ol)计=好(/钞(Kc摸ol绣(j涌)%雷ro建w(锅j)季,j=1,湾nco驶l)芬/)do屑j=懒2,幼nco正lrow韵1=菊lb闭oun膏d(K饥col鹊(j)宏%ro刮w,1欠)!i航_1较jdo雨i伶=衔r更ow伪1,俩j-府1ro胡w_斧1州=夏ma截x(消ro昼w1驴,l姑bo多un弹d(壳Kc古ol臂(i境)%配ro烟w,折1)述)!摄i_阿1s归=奋su欧m(罚di竭ag妨(r性ow豪_1誉:i槐-1痒)战*K摧co振l(园i)似%r序ow互(r滤ow束_1们:i择-1雅)房诚*袖&Kc甲ol纽奉(j堆)%尤ro贡w(脆ro火w_馋1:防i-遥1)演)!求和词部分Kco便l(j汁)%r池ow(抬i)=笑(Kco贤l(j忠)%r抄ow(沉i)-烫s)/紧dia阅g(i)end君do！第j列系数殃分解完低毕s=厦su辞m(d柱iag坝(ro才w1:臭j-1掘)*K歪col决(j)普%ro肉w(r屈ow1部:j-廉1)*血*2)dia楼g(j)=dia蜻g(j)-轮s!第j列的主腾对角元赵素end显do4-5最)一般挥情况