ANSYS 整体刚度矩阵提取子程序——

*deck,rdfull USERDISTRIB ansys
program rdfull
c *** primary function: demonstrates use of binary access routines
c *** secondary function: Read and reformat full file
c
c *** copyright(c) 2009 SAS IP, Inc. All rights reserved.
c *** ansys, inc.
c
c NOTICE - A new assembly process, termed 'symbolic assembly', has 一个被称为symbolic assembly的新的组装过程取代了
c replaced the old assembly process, termed 'frontal 原来被称为 frontal assembly的旧的组装过程，这个
c assembly', and is now the default assembly process for 新的组装过程是现在大多数分析默认的组装过程。
c most analyses. This program demonstrates how to read 本程序展示了怎样去对.FULL文件进行读取和改变格式
c and reformat the .FULL file that was created using 后写入到新的文件，其中这里的.FULL文件可由symbolic
c frontal assembly or symbolic assembly. ANSYS writes the assembly组装过程或frontal assembly组装过程产生。
c .FULL file if the PSOLVE,ELFORM ANSYS将在PSOLVE,ELFORM;PSOLVE,ELPREP;PSOLVE,TRIANG三组命令被依次执行
c PSOLVE,ELPREP 后写出.FULL文件。 当然如果ANSYS使用了sparse求解器、ICCG求解器、
c PSOLVE,TRIANG JCG求解器三者之一或使用了most mode extraction methods也将会
c sequence is used. ANSYS will also write the .FULL file 写出.FULL文件。
c when the sparse, ICCG, or JCG solver is used, as well as
c when most mode extraction methods are used.

c Be sure to set up for modal ANTYPE,2 ？？？？？？
c and Block Lanczos MODOPT,LANB,nmode,0,0, ,OFF ？？？？？？
c (nmode is not used - it can be any value) nmode没有被使用，可以是任何值

c If the free-free stiffness and mass matrices are desired, 如果想获得无约束的刚度矩阵或质量矩阵，则
c make sure there are no constraints on the model. 确保没有对模型施加约束。

#include "impcom.inc"

external ihsort !将数据表中的数据按第某一行中的数据从小到大排列
external binset, binini, binrd8, binclo, biniqr, bintfo
!打开二进制文件，初始化输入输出设备，读取二进制文件中的记录
!关闭打开的二进制文件，二进制文件或系统的查询，定义二进制文件头的必要信息
integer binset, biniqr

integer IOLENG, ROWLENG, WAVEFRONT, MAXNODE, MAXEQN
parameter (IOLENG=16384, ROWLENG=1000, WAVEFRONT=10000,
x MAXNODE=100000, MAXEQN=500000)

LONGINT jloc
integer kunit,munit,units,code,nunit,nbuf,lbuf,npage,keyrw,j,
x kext,i,n,ivect(1

00),kbf,assemb,numdof,lenbac,nontp,
x nmatrx,lumpm,baclst(MAXNODE),nmax,sortlist(2,MAXNODE),
x nmass,nstif,numNodes,nrow,nDofEachNode(MAXEQN),
x dof(MAXEQN),idof,prevDof,currDof,dofMap(MAXEQN),
x irow,node,kdof,lll(ROWLENG),
x mr,nterms,indx(ROWLENG),l(WAVEFRONT)

integer buffer(IOLENG)
double precision krow(ROWLENG), mrow(MAXEQN)

character*80 title(2)
character*106 pname,mname,kname
character*32 jobnam

integer iout,intpdp,lenfnm,reclng

c ***** arrays for reading file data *****
integer iarray(ROWLENG)
double precision darray(ROWLENG)
equivalence (darray(1),iarray(1))



c ********** define the unit numbers **********定义设备编号
iout = 6 !表示输入输出设备为显示器
kunit = 2 !新建文件unit=kunit=2，用来存放整体刚度矩阵信息
munit = 3 !新建文件unit=munit=3，用来存放质量刚度矩阵信息
c ********** define the number of integers per double precision
intpdp = biniqr (0,2) !返回值=2，表示每个双精度数所占的空间等于2个整数所占的空间
c ********** define the number of characters in the file name
lenfnm = biniqr (0,3) !返回值=260，表示所打开的二进制文件的名字最多占用260个字节
c ********** define the i/0 buffer page length (integer*4 words)
reclng = biniqr (0,1)
!返回值为16384，表示标准输入输出缓冲区域页的长度为16384个字符=16kb=系统块的大小
c
c ********** define the file names **********定义文件名
pname = 'file.full' !ANSYS的结果文件，存储了整体刚度矩阵和质量刚度矩阵
mname = 'MASS.MATRIX' !用来存放质量刚度矩阵信息的文件的名字
kname = 'STIFFNESS.MATRIX' !用来存放整体刚度矩阵信息的文件的名字

write (iout,2000)
2000 format (/' ***** WRITE OUT ANSYS MATRICES FROM file.full',
x ' *****'//' MASS MATRIX ON FILE = mass.matrix'/
x ' STIFFNESS MATRIX ON FILE = stiffness.matrix'//
x ' Only the symmetric part of the matrices is written')
!仅仅写出矩阵的对称部分
c
c ********** initialize the bin routines **********
call binini (iout) !初始化输入输出设备
c ********** initialize the header common **********
title(1) = 'New title as given by ROM' !新定义的文件头记录中保存的主标题
title(2) = 'New subtitle as given by ROM' !新定义的文件头记录中保存的副标题
jobnam = 'file ' !新定义的文件头记录中保存的工作名
units = 0 !新定义的文件头记录中保存的单位制
c --- the value for code below should be changed to appropriate 3rd party code
code = 200 !定义

第三方供应商的代码
call bintfo (title,jobnam,units,code)!为子程序binhed8写一个新的文件头记录提供必要的信息
c subroutine bintfo (title,jobnam,units,code)
c *** primary function: set information necessary for binhed
c 为子程序binhed8写一个新的文件头记录提供必要的信息
c --- This routine is intended to be used in standalone programs.
c --- This routine should not be linked into the ANSYS program.
c *** Notice - This file contains ANSYS Confidential information ***
c typ=int,dp,log,chr,dcp siz=sc,ar(n) intent=in,out,inout
c input arguments:
c variable (typ,siz,intent) description
c title (chr*80,ar(2),in) - main title and 1st subtitle 主标题和第一个副标题
c jobnam (chr*8,sc,in) - jobname 文件名
c units (int,sc,in) - units 单位制
c = 0 - user defined units 用户自定义的单位制
c = 1 - SI (MKS)
c = 2 - CSG
c = 3 - U.S. Customary, using feet
c = 4 - U.S. Customary, using inches
c = 6 - MPA
c = 7 - uMKS
c code (int,sc,in) - code defining 3rd party vendor 定义第三方供应商的代码
c (contact ANSYS, Inc. for code assignment)
c
c output arguments:
c none
c
c ********** open the ASCII files **********新建有格式顺序文件
open (unit=kunit,file=kname,status='unknown')
!新建文件unit=kunit=2，file=kname='STIFFNESS.MATRIX'，用来存放整体刚度矩阵信息
rewind 2 !将文件2记录的读写位置定位在文件的起点
open (unit=munit,file=mname,status='unknown')
!新建文件unit=munit=3，file=mname='MASS.MATRIX'，用来存放质量刚度矩阵信息
rewind 3 !将文件3记录的读写位置定位在文件的起点
c
c ********** open the file **********
c
c ***** set the i/o unit number
nunit =7 !为打开的二进制文件file.full设置一个编号
c ***** define buffer number *****
nbuf = 1 !为打开的二进制文件所存放的缓冲区域设置一个编号
c ***** start address in the buffer *****
lbuf = 1 !用工作数组的列号来表示的缓冲区域开始的位置
c ***** number of pages for nbuf *****
npage = 1 !设置编号为nbuf的缓冲区域所包含的缓冲区域页的数量
c ***** read key *****
keyrw = 1 !文件读写权限中，读的编号
c ***** length of each page *****
j = reclng !缓冲区域页的长度=系统记录的长度=系统块的大小=16384字节
c ***** read external file format *****
kext = 1 !读取外部文件的格式
c
i = binset (nbuf,nunit,keyrw,lbuf,j,npage,pname,lenfnm,kext,
x buffer(1))
c function binset (nblk,nunit,ikeyrw,istart,paglen,npage,pname,nchar,kext,Buffer4)
c *** primary function: initialize paging s

pace for a blocked binary file.
c 为一个被分块化的二进制文件分配缓冲区域即对即将存放二进制文件的
c 缓冲区域（等于所有缓冲区域页的和）进行初始化
c binset should be used to open a blocked file
c binset程序的功能是打开一个被分块化的二进制文件
c before binrd8 or binwrt8 are used. binclo should
c be used to close the file.
c --- This routine is intended to be used in standalone programs.
c --- This routine should not be linked into the ANSYS program.
c *** Notice - This file contains ANSYS Confidential information ***
c input arguments:
c nblk (int,sc,in) - block number (1 to BIO_MAXFILES max)
c 二进制文件所存放的缓冲区域的编号
c 一个缓冲区域=其下所有缓冲区域页的和
c nunit (int,sc,in) - fortran unit number for the file
c (if 0, bit bucket)
c 打开的二进制文件的编号
c ikeyrw (int,sc,in) - read/write flag 设置读写权限
c = 0 - both read & write
c = 1 - read
c = 2 - write
c = 9 - read only
c istart (int,sc,in) - starting location in buffer array
c 用工作数组的列号来表示的缓冲区域开始的位置
c usually 1 for nblk=1,
c paglen*npage+1 for nblk=2,etc.
c paglen (int,sc,in) - page length in integer*4 words for external
c files 缓冲区域页的长度=系统记录的长度=系统块的大小=16384字节
c paglen should always be a multiple of
c 512 words for efficiency
c npage (int,sc,in) - number of pages (1 to BIO_MAXBLOCKS max) 缓冲区域所包含的缓冲区域页的数量
c pname (chr,ar(*),in) - name of the file 打开的二进制文件的名字
c nchar (int,sc,in) - number of characters in the file name(not used)
c 表示所打开的二进制文件的名字最多占用的字节数量
c kext (int,sc,in) - no longer used, always external format
c 文件的扩展名的编号
c Buffer4 (i4, ar(*),inout) - work array for paging, should be
c dimensioned to paglen*npage*nblk (max)
c 读取二进制文件时临时使用的数组，应该被定义为paglen*npage*nblk (max)维数
c output arguments:
c binset (int,func,out) - error status 文件打开的状态
c = 0 - no error 返回值等于0，则文件顺利打开
c <>0 - error occurred 返回值不等于0，则文件打开失败
c Buffer4 (i4, ar(*),inout) - work array for paging
c 读取二进制文件时临时使用的数组
c
c ***** check for error *****
if (i/=0) then
write

(iout,2004) pname
2004 format (/' *** ERROR ***'/' DATA FILE ',a,' DOES NOT EXIST ')
go to 999
endif

c ***** file header *****
n = 100
jloc = 0
call binrd8 (nbuf,jloc,n,ivect(1),kbf,buffer(1)) !读取二进制文件的第一条记录
!从缓冲区域读取二进制文件，一次只读取一条记录
c subroutine binrd8 (nblk,LongLocL,leng,ivect,kbfint,Buffer4)
c ********** buffer read routine **********
c *** Notice - This file contains ANSYS Confidential information ***
c input arguments:
c nblk (int,sc,in) - block number. see fd___(i.e. fdtri for tri
c 二进制文件所存放的缓冲区域的编号,应与binset程序中的编号一致
c 一个缓冲区域=其下所有缓冲区域页的和
c (as defined with subroutine bioset)
c LongLocL(LONG,sc,inout)- location in integer*4 words of the startin
c position on the file.
c 在被分块的二进制文件中当前读取数据的位置
c leng (int,sc,inout) - number of words to read into ivect. (must be
c less or equal to dimension given to ivect in
c the calling routine). if ivect is to be used
c as integers, use as is. if ivect is to be
c used for double precision numbers, it must be
c increased by multiplying it by INTPDP（一个全局变量）.
c if negative, skip record and do not return
c data(results).
c data(results).
c 为即将读取的数据的个数设置一个初始值<=buffer数组的维数,这个参数不能忽略
c Buffer4 (i4, ar(*),inout) - work array for paging, should be the
c same array as used in binset
c 读取二进制文件时临时使用的数组
c output arguments:
c LongLocL(LONG,sc,inout)- location in integer*4 words of the current
c position on the file. It is updated after
c each read operation
c 读取完一条记录后即将读取下一条记录的数据的位置，以4个字节为一个单位
c leng (int,sc,inout) - tells you how many items it actually read(in
c integer words).一次所实际读取到的数据的数量，以4个字节为
c 一个单位，如果没有读到数据则n=0，故一个双精度数据算2个数据
c if zero, end of file(error case)
c ivect (int,ar(*),out) - results (can be either integer or double
c precision in the calling routine)
c 实际所读取到的数据（可以是整型或双精度型）依次存储在这个数组里面
c kbfint (int,sc,out) - key for type(used only for AUX2 dump) 读取到的数据类型
c = 0 double preci

sion data 读取到的数据类型为双精度
c > 0 integer data(usually the same as leng)
c 读取到的数据类型为整型，并返回读取到的整型数据的数量
c Buffer4 (i4,ar(*),inout) - work array for paging 读取二进制文件时临时使用的数组
n = 40
call binrd8 (nbuf,jloc,n,ivect(1),kbf,buffer(1)) !读取二进制文件的第二条记录

assemb = ivect(1) !文件的编号
numdof = ivect(8) !每个结点的自由度数量
lenbac = ivect(7) !结点的数量
nontp = ivect(2) !文件中所含方程式的数量
nmatrx = ivect(4) !文件中所含矩阵的数量
lumpm = ivect(11) !质量块的标识
write (iout,2001) numdof,lenbac,lumpm
2001 format(/' Number of DOF per node=',i6/
x ' Number of nodes =',i6/
x ' Lumped mass flag =',i6)

n = numdof !每个结点的自由度数量
call binrd8 (nbuf,jloc,n,ivect(1),kbf,buffer(1)) !读取二进制文件的第三条记录

n = lenbac !结点的数量
call binrd8 (nbuf,jloc,n,baclst(1),kbf,buffer(1))
!读取二进制文件的第四条记录,第四条记录依次存储的是结点的编号，如3，8，2，9，22
!可以用write(*,"(6i6)")(baclst(i),i=1,lenbac)语句来查看结点号的实际排列顺序
!一般都不是按大小顺序排列的

c ********** compress nodes into grid point order *********
c 将所有的结点号压缩并按原来的次序输出，如3，8，2，9，22经过压缩后成为2，3，1，4，5
nmax = 0
do i = 1,lenbac
sortlist(1,i) = baclst(i) !使数组sortlist的第一列依次存储的是结点的编号
sortlist(2,i) = i !使数组sortlist的第二列依次存储的是1,2,3....lenbac(结点的数量)
if (baclst(i)>nmax) nmax = baclst(i) !nmax为结点的最大编号
enddo
write (iout,2002) nmax !结点编号的最大值
2002 format (' Maximum node number =',i6)
call ihsort (2,lenbac,sortlist(1,1),1)
!将数组sortlist中的数据按第1行中的数据从小到大排列
c subroutine ihsort (nper,n,table,iloc) 将表table中的数据按第iloc行中的数据从小到大排列
c *** primary function: quick heap sort of an integer table
c input arguments:
c nper (int,sc,in) - number of terms 数据表table的总行数
c n (int,sc,in) - number of items in the 数据表table中的总列数
c table (int,ar(nper,n),inout) - the table to be sorted 需要整理的数据表
c iloc (int,sc,in) - location in table for sort 按第iloc行中的数据从小到大排列
c output arguments:
c table (int,ar(nper,n),inout) - the sorted table 返回整理好的数据表
do i = 1,lenbac
sortlist(1,i) = i !使数组sortlist的第一列依次存储的是1,2,3....l

enbac(结点的数量)
enddo
call ihsort (2,lenbac,sortlist(1,1),2)
!将数组sortlist中的数据按第2行中的数据从小到大排列
do i = 1,lenbac
baclst(i) = sortlist(1,i) !使数组baclst依次存放的是经过压缩后的结点编号
enddo

c ***** put headers on the files *****
write (kunit,2010)
2010 format ('# STIFFNESS.MATRIX'/'#'/'# FROM ANSYS .FULL FILE'/'#')
write (munit,2011)
2011 format ('# MASS.MATRIX'/'#'/'# FROM ANSYS .FULL FILE'/'#')

c ***** jump to proper code depeding on assembly process *****!判断矩阵的构建方式
if (assemb==4) goto 300 !矩阵按frontal assembly方式建立
if (assemb==-4) goto 600 !矩阵按symbolic assembly方式建立
goto 999

c ***** code to read and write frontal assembly .FULL file *****
300 nmass = 0
nstif = 0
nrow = 0
do irow = 1,nontp

n = 10
call binrd8 (nbuf,jloc,n,lll(1),kbf,buffer(1))
mr = lll(1)
kdof = abs(lll(2))
node = (kdof-1)/numdof
idof = kdof - node*numdof
node = baclst(node+1)
kdof = (node-1)*numdof + idof
do i = 3,n
nrow = nrow + 1
l(nrow) = lll(i)
enddo

nterms = 10
call binrd8 (nbuf,jloc,nterms,indx(1),kbf,buffer(1))
do i = 1,nterms
idof = l(indx(i))
node = (idof-1)/numdof
idof = idof - node*numdof
node = baclst(node+1)
lll(i) = (node-1)*numdof + idof
enddo
n = 10
call binrd8 (nbuf,jloc,n,krow(1),kbf,buffer(1))
i = 1
if (lumpm .ge. 1) i = 2
n = n/intpdp-i
do i = 1,n
if (krow(i) .ne. 0.0d0) then
nstif = nstif + 1
if (kdof .lt. lll(i)) then
write (kunit,2005) kdof,lll(i),krow(i)
else
write (kunit,2005) lll(i),kdof,krow(i)
endif
endif
enddo

if (nmatrx.gt.1 .and. lumpm.eq.0) then
n = 10
call binrd8 (nbuf,jloc,n,mrow(1),kbf,buffer(1))
n = n/intpdp
do i = 1,n
if (mrow(i) .ne. 0.0d0) then
nmass = nmass + 1
if (kdof .lt. lll(i)) then
write (munit,2005) kdof,lll(i),mrow(i)
2005 format (2i8,1pe20.12)
else
write (munit,2005) lll(i),kdof,mrow(i)
endif
endif
enddo
else
mrow(1) = krow(n+2)
if (mrow(1) .ne. 0.0d0) then
nmass = nmass + 1
write (munit,2005) kdof,kdof,mrow(1)
endif
endif

l(mr) = l(nrow)
l(nrow) = 0
nrow = nrow - 1

enddo

clo

se (unit=kunit,status='keep')
close (unit=munit,status='keep')

write (iout,2008) nmass, nstif
2008 format (/' Number of mass matrix terms =',i8/
x ' Number of stiffness matrix terms=',i8)
goto 999


c ***** code to read and write symbolic assembled .FULL file *****
600 nmass = 0
nstif = 0
nrow = 0
numNodes = ivect(33) !被包含在矩阵内部的结点的数量

c *** read nodal extent vector 读取每个结点的自由度数量
n = lenbac !结点的数量
jloc = ivect(36) !指向自由度信息的指针
call binrd8 (nbuf,jloc,n,nDofEachNode(1),kbf,buffer(1))
!读取结点自由度信息的记录，返回各个结点的自由度数量
!与数组baclst依次存放的经过压缩后的结点编号相对应
c *** read dof vector
n = nontp !文件中所含方程式的数量
call binrd8 (nbuf,jloc,n,dof(1),kbf,buffer(1))
!读取自由度维数信息的记录，返回各个结点的自由度的编号
!如1,2,3,4,5,6等，如果返回的值为负值，则代表这个自由度被约束住
!与数组baclst依次存放的经过压缩后的结点编号相对应
!返回的n的值等于每个结点的自由度的和
!数组dof中返回的是矩阵中某一行所对应的结点自由度的编号
idof = 0
prevDof = 0
do i = 1,numNodes !被包含在矩阵内部的结点的数量
!i代表第几个结点
do j = 1,nDofEachNode(i) !nDofEachNode(i)代表第i个结点的自由度数量
idof = idof + 1
currDof = dof(idof) !currDof代表第i个结点的第j个自由度的编号
if (currDof>0) then
dofMap(idof) = prevDof + currDof
!dofMap返回的是没有被约束住的i结点的第j个自由度在矩阵中对应的行号
else
dofMap(idof) = -(prevDof + abs(currDof))
endif
enddo
prevDof = prevDof + numdof !numdof为每个结点的自由度数量
enddo


!读取整体刚度矩阵,

jloc = ivect(19) !指向整体刚度矩阵的指针
do irow = 1,nontp !文件中所含方程式的数量
n = 10
call binrd8 (nbuf,jloc,n,lll(1),kbf,buffer(1))
write(*,"(6i6)")(lll(i),i=1,n)
!返回矩阵第irow行中数据/=0的列号
!***************其中最后一个是矩阵的irow行对角线上的列号
i = n*intpdp
call binrd8 (nbuf,jloc,i,krow(1),kbf,buffer(1))
!返回矩阵第irow行中数据/=0的列号所对应的具体数据（双精度）
if (dof(irow)<=0) cycle !数组dof中返回的是矩阵中某一行所对应的结点自由度的编号
!不读取某结点的某个被约束住的自由度所对应的矩阵行，因为这一行所表示的方程对求解没有用处
idof = dofMap(irow) !实际上dofMap(idof)=矩

阵的行号
node = (idof-1)/numdof
idof = idof - node*numdof !返回的是第irow行所对应的结点自由度编号
node = baclst(node+1) !返回的是第irow行所对应的结点的实际编号
kdof = (node-1)*numdof + idof !如果刚度矩阵对应的结点是按节点号从小到大的顺序排列的
!******************则kdof代表以上读取到的数据在这种刚度矩阵中所处的是第kdof行
!注意之前所讨论的刚度矩阵都是与数组baclst依次存放的经过压缩后的结点编号相对应
!而实际的刚度矩阵对应的结点是按节点号从小到大的顺序排列的，且具有对称性
do i = 1,n
idof = dofMap(lll(i)) !列号lll(i)所对应的矩阵的行号
node = (idof-1)/numdof
idof = idof - node*numdof !返回的是第lll(i)列所对应的结点自由度编号
node = baclst(node+1) !返回的是第lll(i)列所对应的结点的实际编号
lll(i) = (node-1)*numdof + idof !如果刚度矩阵对应的结点是按结点号从小到大的顺序排列的
!***************则lll(i)代表以上读取到的数据在这种刚度矩阵中所处的是第lll(i)列
!注意之前所讨论的刚度矩阵都是与数组baclst依次存放的经过压缩后的结点编号相对应
!而实际的刚度矩阵对应的结点是按节点号从小到大的顺序排列的，且具有对称性
if (krow(i)/=0.0d0) then
nstif = nstif + 1
if (kdof!.FULL文件中保存的数据量仅是整体刚度矩阵的一半的数据量即上三角或下三角的数据量，
!而且这些数据都在与数组baclst依次存放的经过压缩后的结点编号相对的刚度矩阵的上三角的位置上。
!因此如果要得到与按节点号从小到大的顺序排列的结点相对应的刚度矩阵，
!需要采用一定的方法使这些数据都在这样的矩阵的上三角或下三角被输出。
write (kunit,2105) kdof,lll(i),krow(i)
else !kdof,lll(i)分别表示数据krow(i)所在的行号和列号
write (kunit,2105) lll(i),kdof,krow(i)
endif
endif
enddo
enddo


!读取质量刚度矩阵

if (nmatrx>1) then
jloc = ivect(27) !指向质量刚度矩阵的指针
if (lumpm .eq. 0) then
do irow = 1,nontp
n = 10
call binrd8 (nbuf,jloc,n,lll(1),kbf,buffer(1))
i = n*intpdp
call binrd8 (nbuf,jloc,i,mrow(1),kbf,buffer(1))
if (dof(irow) .le. 0) cycle
idof = dofMap(irow)
node = (idof-1)/numdof
idof = idof - node*numdof
node = baclst(node+1)

kdof = (node-1)*numdof + idof
do i = 1,n
if (dof(lll(i)) .le. 0) cycle
idof = dofMap(lll(i))
node = (idof-1)/numdof
idof = idof - node*numdof
node = baclst(node+1)
lll(i) = (node-1)*numdof + idof
if (mrow(i) .ne. 0.0d0) then
nmass = nmass + 1
if (kdof .lt. lll(i)) then
write (munit,2105) kdof,lll(i),mrow(i)
2105 format (2i8,1pe20.12)
else
write (munit,2105) lll(i),kdof,mrow(i)
endif
endif
enddo
enddo
else
i = nontp*intpdp
call binrd8 (nbuf,jloc,i,mrow(1),kbf,buffer(1))
do irow = 1,nontp
if (dof(irow) .le. 0) cycle
idof = dofMap(irow)
node = (idof-1)/numdof
idof = idof - node*numdof
node = baclst(node+1)
kdof = (node-1)*numdof + idof
if (mrow(irow) .ne. 0.0d0) then
nmass = nmass + 1
write (munit,2105) kdof,kdof,mrow(irow)
endif
enddo
endif
endif

!关闭并保存用来存储刚度矩阵和质量矩阵信息的两个有格式顺序文件
close (unit=kunit,status='keep')
close (unit=munit,status='keep')

write (iout,2108) nmass, nstif
2108 format (/' Number of mass matrix terms =',i8/
x ' Number of stiffness matrix terms=',i8)

!关闭并保存二进制文件.FULL
999 call binclo (nbuf,'KEEP',buffer(1))
end


!子程序ihsort的建立

*deck,ihsort parallel
subroutine ihsort (nper,n,table,iloc)
c *** primary function: quick heap sort of an integer table

c input arguments:
c nper (int,sc,in) - number of terms
c n (int,sc,in) - number of items in the table
c table (int,ar(nper,n),inout) - the table to be sorted
c iloc (int,sc,in) - location in table for sort

c output arguments:
c table (int,ar(nper,n),inout) - the sorted table


integer nper, n, table(nper,n), iloc, i,j,k,l,m, hold, ip

m = n
10 m = m/2
if (m==0) go to 50
k = n - m
j = 1
20 i = j
30 l = i + m
if (table(iloc,i)<=table(iloc,l)) go to 40
do ip = 1,nper
hold = table(ip,i)
table(ip,i) = table(ip,l)
table(ip,l) = hold
enddo
i = i - m
if (i>=1) go to 30
40 j = j + 1
if (j>k) go to 10
go to 20

50 return
end


c
c 提取整体刚度矩阵程

序的基本结构：
c 变量声明
c 初始化输入输出设备
c 为写二进制文件头记录定义必要的信息
c 新建两个有格式顺序文件用于存储整体刚度矩阵和质量矩阵
c 打开二进制文件.FULL
c 对读取二进制文件中的相关数据写出到新建的两个有格式顺序文件中
c 关闭并保存两个有格式顺序文件
c 关闭并保存已打开的二进制文件

——Jake-Ni