Context Navigation

source: issm/trunk-jpl/src/c/modules/AllocateSystemMatricesx/AllocateSystemMatricesx.cpp@ 23641

Last change on this file since 23641 was 23641, checked in by Mathieu Morlighem, 6 years ago
CHG: speeding up allocation without MPI_Allreduce
File size: 9.0 KB

Rev	Line
[16126]	1	/*!\file AllocateSystemMatricesx
	2	* \brief retrieve vector from inputs in elements
	3	*/
	4
	5	#include "./AllocateSystemMatricesx.h"
	6	#include "../../shared/shared.h"
	7	#include "../../toolkits/toolkits.h"
	8
	9	void AllocateSystemMatricesx(Matrix<IssmDouble> pKff,Matrix<IssmDouble> pKfs,Vector<IssmDouble> pdf,Vector<IssmDouble> ppf,FemModel* femmodel){
	10
	11	/Intermediary/
	12	int fsize,ssize,flocalsize,slocalsize;
	13	int connectivity, numberofdofspernode;
	14	int m,n,M,N;
	15	int *d_nnz = NULL;
	16	int *o_nnz = NULL;
	17
	18	/output/
	19	Matrix<IssmDouble> *Kff = NULL;
	20	Matrix<IssmDouble> *Kfs = NULL;
	21	Vector<IssmDouble> *pf = NULL;
	22	Vector<IssmDouble> *df = NULL;
	23
	24	bool oldalloc=false;
[16438]	25	char* toolkittype=NULL;
[16126]	26
	27	/retrieve parameters: /
	28	femmodel->parameters->FindParam(&connectivity,MeshAverageVertexConnectivityEnum);
	29
	30	/retrieve node info/
[23587]	31	fsize = femmodel->nodes->NumberOfDofs(FsetEnum);
	32	ssize = femmodel->nodes->NumberOfDofs(SsetEnum);
	33	flocalsize = femmodel->nodes->NumberOfDofsLocal(FsetEnum);
	34	slocalsize = femmodel->nodes->NumberOfDofsLocal(SsetEnum);
[16126]	35
[23587]	36	numberofdofspernode=femmodel->nodes->MaxNumDofs(GsetEnum);
[16126]	37
[16438]	38	/*if our matrices are coming from issm, we don't do dynamic allocation like Petsc
[16441]	39	* does, and this routine is essentially useless. Force standard alloc in this case: */
[16438]	40	toolkittype=ToolkitOptions::GetToolkitType();
	41
[16126]	42	if(oldalloc){
	43	if(pKff) Kff=new Matrix<IssmDouble>(fsize,fsize,connectivity,numberofdofspernode);
	44	if(pKfs) Kfs=new Matrix<IssmDouble>(fsize,ssize,connectivity,numberofdofspernode);
	45	if(pdf) df =new Vector<IssmDouble>(fsize);
	46	if(ppf) pf =new Vector<IssmDouble>(fsize);
	47	}
	48	else{
	49	if(pKff){
	50	m=flocalsize; n=flocalsize; /local sizes/
	51	M=fsize; N=fsize; /global sizes/
[16441]	52	if(strcmp(toolkittype,"issm")==0){
	53	Kff=new Matrix<IssmDouble>(m,n,M,N,NULL,NULL);
	54	}
	55	else{
	56	MatrixNonzeros(&d_nnz,&o_nnz,femmodel,FsetEnum,FsetEnum);
	57	Kff=new Matrix<IssmDouble>(m,n,M,N,d_nnz,o_nnz);
	58	xDelete<int>(d_nnz);
	59	xDelete<int>(o_nnz);
	60	}
[16126]	61	}
	62	if(pKfs){
	63	m=flocalsize; n=slocalsize; /local sizes/
	64	M=fsize; N=ssize; /global sizes/
[16441]	65	if(strcmp(toolkittype,"issm")==0){
	66	Kfs=new Matrix<IssmDouble>(m,n,M,N,NULL,NULL);
	67	}
	68	else{
	69	MatrixNonzeros(&d_nnz,&o_nnz,femmodel,FsetEnum,SsetEnum);
	70	Kfs=new Matrix<IssmDouble>(m,n,M,N,d_nnz,o_nnz);
	71	xDelete<int>(d_nnz);
	72	xDelete<int>(o_nnz);
	73	}
[16126]	74	}
	75	if(pdf) df =new Vector<IssmDouble>(flocalsize,fsize);
	76	if(ppf) pf =new Vector<IssmDouble>(flocalsize,fsize);
	77	}
[23066]	78
[16438]	79	/Free ressources: /
	80	xDelete<char>(toolkittype);
[16126]	81
	82	/Allocate output pointers/
	83	if(pKff) *pKff = Kff;
	84	if(pKfs) *pKfs = Kfs;
	85	if(pdf) *pdf = df;
	86	if(ppf) *ppf = pf;
	87	}
	88
	89	void MatrixNonzeros(int pd_nnz,int po_nnz,FemModel* femmodel,int set1enum,int set2enum){
	90
	91	/Intermediary/
	92	int i,j,k,index,offset,count;
	93	int d_nz,o_nz;
	94	Element *element = NULL;
	95	Load *load = NULL;
	96	int *head_e = NULL;
	97	int *next_e = NULL;
	98	int *count2offset_e = NULL;
	99	int *head_l = NULL;
	100	int *next_l = NULL;
	101	int *count2offset_l = NULL;
	102	int *lidlist = NULL;
	103
	104	/output/
	105	int *d_nnz = NULL;
	106	int *o_nnz = NULL;
	107
	108	/Get vector size and number of nodes/
[23587]	109	int numnodes = femmodel->nodes->NumberOfNodes();
[23641]	110	int localmasters = femmodel->nodes->NumberOfNodesLocal();
[16126]	111	int localnumnodes = femmodel->nodes->Size();
[23587]	112	int numberofdofspernode = femmodel->nodes->MaxNumDofs(GsetEnum);
	113	int M = femmodel->nodes->NumberOfDofs(set1enum);
	114	int N = femmodel->nodes->NumberOfDofs(set2enum);
	115	int m = femmodel->nodes->NumberOfDofsLocal(set1enum);
	116	int n = femmodel->nodes->NumberOfDofsLocal(set2enum);
[16126]	117	int numnodesperelement = femmodel->elements->MaxNumNodes();
[23587]	118	int numnodesperload = femmodel->loads->MaxNumNodes();
[16126]	119
	120	/*First, we are building chaining vectors so that we know what nodes are
	121	* connected to what elements. These vectors are such that:
	122	* for(int i=head[id];i!=-1;i=next[i])
	123	* will loop over all the elements that are connected to the node number
	124	* id*/
	125	head_e = xNew<int>(localnumnodes); for(i=0;i<localnumnodes;i++) head_e[i]=-1;
	126	next_e = xNew<int>(femmodel->elements->Size()*numnodesperelement);
	127	count2offset_e = xNew<int>(femmodel->elements->Size()*numnodesperelement);
	128
	129	k=0;
	130	for(i=0;i<femmodel->elements->Size();i++){
[18521]	131	element = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
[16126]	132	lidlist = xNew<int>(element->GetNumberOfNodes());
	133	element->GetNodesLidList(lidlist);
	134
	135	for(j=0;j<element->GetNumberOfNodes();j++){
	136	index = lidlist[j];
	137	_assert_(index>=0 && index<numnodes);
	138
	139	count2offset_e[k]=i;
	140	next_e[k]=head_e[index];
	141	head_e[index]=k++;
	142	}
	143	for(j=0;j<numnodesperelement-element->GetNumberOfNodes();j++) k++;
	144
	145	xDelete<int>(lidlist);
	146	}
	147
	148	/Chain for loads/
	149	head_l = xNew<int>(localnumnodes); for(i=0;i<localnumnodes;i++) head_l[i]=-1;
[23587]	150	next_l = xNew<int>(femmodel->loads->Size()*numnodesperload);
	151	count2offset_l = xNew<int>(femmodel->loads->Size()*numnodesperload);
[16126]	152	k=0;
	153	for(i=0;i<femmodel->loads->Size();i++){
[18521]	154	load = xDynamicCast<Load*>(femmodel->loads->GetObjectByOffset(i));
[16126]	155	lidlist = xNew<int>(load->GetNumberOfNodes());
	156	load->GetNodesLidList(lidlist);
	157
	158	for(j=0;j<load->GetNumberOfNodes();j++){
	159	index = lidlist[j];
	160	_assert_(index>=0 && index<numnodes);
	161
	162	count2offset_l[k]=i;
	163	next_l[k]=head_l[index];
	164	head_l[index]=k++;
	165	}
	166	for(j=0;j<numnodesperload-load->GetNumberOfNodes();j++) k++;
	167
	168	xDelete<int>(lidlist);
	169	}
	170
	171	/OK now count number of dofs and flag each nodes for each node i/
	172	bool *flags = xNew<bool>(localnumnodes);
	173	int *flagsindices = xNew<int>(localnumnodes);
[23641]	174	int *d_connectivity = xNewZeroInit<int>(localnumnodes);
	175	int *o_connectivity = xNewZeroInit<int>(localnumnodes);
[16126]	176	int *all_connectivity_clone = xNewZeroInit<int>(numnodes);
	177
[23641]	178	Vector<IssmDouble> *connectivity_clone= new Vector<IssmDouble>(localmasters,numnodes);
	179
[16126]	180	/*Resetting flags to false at eahc iteration takes a lot of time, so we keep track of the flags
	181	* to reset in flagsindices, initialized with -1*/
	182	for(i = 0;i<localnumnodes;i++) flags[i] = false;
	183	for(i = 0;i<localnumnodes;i++) flagsindices[i] = -1;
	184
	185	/Create connectivity vector/
	186	for(i=0;i<femmodel->nodes->Size();i++){
[18521]	187	Node* node=xDynamicCast<Node*>(femmodel->nodes->GetObjectByOffset(i));
[16126]	188
[23499]	189	/Reinitialize flags to false/
	190	j=0;
	191	while(j<localnumnodes){
	192	if(flagsindices[j]>=0){
	193	flags[flagsindices[j]] = false;
	194	flagsindices[j] = -1;
	195	j++;
[16126]	196	}
[23499]	197	else{
	198	break;
	199	}
	200	}
[16126]	201
[23499]	202	//for(j=0;j<localnumnodes;j++) flags[j]=false;
[16126]	203
[23499]	204	/Loop over elements that hold node number i/
	205	//if(head_e[node->Lid()]==-1 && head_l[node->Lid()]==-1){
	206	// printf("[%i] vertex %i\n",IssmComm::GetRank(),node->Lid()+1);
	207	//}
	208	for(j=head_e[node->Lid()];j!=-1;j=next_e[j]){
	209	offset=count2offset_e[j];
	210	element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(offset));
	211	element->SetwiseNodeConnectivity(&d_nz,&o_nz,node,flags,flagsindices,set1enum,set2enum);
	212	if(node->IsClone()){
[23641]	213	connectivity_clone->SetValue(node->Pid(),d_nz+o_nz,ADD_VAL);
[16126]	214	}
[23499]	215	else{
[23641]	216	d_connectivity[node->Lid()]+=d_nz;
	217	o_connectivity[node->Lid()]+=o_nz;
[16126]	218	}
	219	}
[23499]	220	for(j=head_l[node->Lid()];j!=-1;j=next_l[j]){
	221	offset=count2offset_l[j];
	222	load=xDynamicCast<Load*>(femmodel->loads->GetObjectByOffset(offset));
	223	load->SetwiseNodeConnectivity(&d_nz,&o_nz,node,flags,flagsindices,set1enum,set2enum);
	224	if(node->IsClone()){
[23641]	225	connectivity_clone->SetValue(node->Pid(),d_nz+o_nz,ADD_VAL);
[23499]	226	}
	227	else{
[23641]	228	d_connectivity[node->Lid()]+=d_nz;
	229	o_connectivity[node->Lid()]+=o_nz;
[23499]	230	}
	231	}
[16126]	232	}
	233	xDelete<bool>(flags);
	234	xDelete<int>(flagsindices);
	235	xDelete<int>(count2offset_e);
	236	xDelete<int>(head_e);
	237	xDelete<int>(next_e);
	238	xDelete<int>(count2offset_l);
	239	xDelete<int>(head_l);
	240	xDelete<int>(next_l);
	241
	242	/sum over all cpus/
[23641]	243	connectivity_clone->Assemble();
	244	IssmDouble* serial_connectivity_clone=NULL;
	245	femmodel->GetLocalVectorWithClonesVertices(&serial_connectivity_clone,connectivity_clone);
	246	delete connectivity_clone;
[16126]	247
	248	if(set1enum==FsetEnum){
	249	count=0;
	250	d_nnz=xNew<int>(m);
	251	o_nnz=xNew<int>(m);
	252	for(i=0;i<femmodel->nodes->Size();i++){
[18521]	253	Node* node=xDynamicCast<Node*>(femmodel->nodes->GetObjectByOffset(i));
[23499]	254	if(!node->IsClone()){
[23612]	255	for(j=0;j<node->fsize;j++){
[16126]	256	_assert_(count<m);
[23641]	257	d_nnz[count]=numberofdofspernode*(d_connectivity[node->Lid()] + reCast<int>(serial_connectivity_clone[node->Lid()]));
	258	o_nnz[count]=numberofdofspernode*(o_connectivity[node->Lid()] + reCast<int>(serial_connectivity_clone[node->Lid()]));
[16126]	259	if(d_nnz[count]>n) d_nnz[count]=n;
	260	if(o_nnz[count]>N-n) o_nnz[count]=N-n;
	261	count++;
	262	}
	263	}
	264	}
	265	_assert_(m==count);
	266	}
	267	else{
	268	_error_("STOP not implemented");
	269	}
	270	xDelete<int>(d_connectivity);
	271	xDelete<int>(o_connectivity);
[23641]	272	xDelete<IssmDouble>(serial_connectivity_clone);
[16126]	273
	274	/Allocate ouptput pointer/
	275	*pd_nnz=d_nnz;
	276	*po_nnz=o_nnz;
	277	}

Note: See TracBrowser for help on using the repository browser.

Download in other formats: