1 /*
2  *  This is an easy to call version of the sdp routine.  It takes as
3  *  input a problem (n,k,C,a,constraints,constant_offset), and an
4  *  initial solution (X,y,Z), allocates working storage, and calls sdp()
5  *  to solve the problem.  The solution is returned in X,y,Z,pobj,dobj, and
6  *  the return code from sdp is returned as the return value from easy_sdp.
7  *
8  */
9 
10 #include <stdlib.h>
11 #include <stdio.h>
12 #include <math.h>
13 #include "declarations.h"
14 
easy_sdp(n,k,C,a,constraints,constant_offset,pX,py,pZ,ppobj,pdobj)15 int easy_sdp(n,k,C,a,constraints,constant_offset,pX,py,pZ,ppobj,pdobj)
16      int n;
17      int k;
18      struct blockmatrix C;
19      double *a;
20      struct constraintmatrix *constraints;
21      double constant_offset;
22      struct blockmatrix *pX;
23      double **py;
24      struct blockmatrix *pZ;
25      double *ppobj;
26      double *pdobj;
27 {
28   int ret;
29   struct constraintmatrix fill;
30   struct paramstruc params;
31   struct blockmatrix work1;
32   struct blockmatrix work2;
33   struct blockmatrix work3;
34   struct blockmatrix bestx;
35   struct blockmatrix bestz;
36   struct blockmatrix Zi;
37   struct blockmatrix dZ;
38   struct blockmatrix dX;
39   struct blockmatrix cholxinv;
40   struct blockmatrix cholzinv;
41   double *workvec1;
42   double *workvec2;
43   double *workvec3;
44   double *workvec4;
45   double *workvec5;
46   double *workvec6;
47   double *workvec7;
48   double *workvec8;
49   double *diagO;
50   double *Fp;
51   double *O;
52   double *dy;
53   double *dy1;
54   double *rhs;
55   double *besty;
56   int printlevel;
57   int ldam;
58   struct sparseblock **byblocks;
59   struct sparseblock *ptr;
60   struct sparseblock *oldptr;
61   int i;
62   int j;
63   struct sparseblock *p;
64   struct sparseblock *q;
65   struct sparseblock *prev=NULL;
66   double gap;
67   int nnz;
68   int denseblocks;
69   int numblocks;
70 
71    /*
72     *  Initialize the parameters.
73     */
74 
75    initparams(&params,&printlevel);
76 
77 
78   /*
79    *  Allocate working storage
80    */
81 
82   alloc_mat(C,&work1);
83   alloc_mat(C,&work2);
84   alloc_mat(C,&work3);
85   alloc_mat_packed(C,&bestx);
86   alloc_mat_packed(C,&bestz);
87   alloc_mat_packed(C,&cholxinv);
88   alloc_mat_packed(C,&cholzinv);
89 
90   besty=(double *)malloc(sizeof(double)*(k+1));
91    if (besty == NULL)
92      {
93        printf("Storage Allocation Failed!\n");
94        exit(205);
95      };
96 
97    if (n > k)
98      {
99        workvec1=(double *)malloc(sizeof(double)*(n+1));
100      }
101    else
102      {
103        workvec1=(double *)malloc(sizeof(double)*(k+1));
104      };
105    if (workvec1 == NULL)
106      {
107        printf("Storage Allocation Failed!\n");
108        exit(205);
109      };
110 
111    if (n > k)
112      {
113        workvec2=(double *)malloc(sizeof(double)*(n+1));
114      }
115    else
116      {
117        workvec2=(double *)malloc(sizeof(double)*(k+1));
118      };
119    if (workvec2 == NULL)
120      {
121        printf("Storage Allocation Failed!\n");
122        exit(205);
123      };
124 
125    if (n > k)
126      {
127        workvec3=(double *)malloc(sizeof(double)*(n+1));
128      }
129    else
130      {
131        workvec3=(double *)malloc(sizeof(double)*(k+1));
132      };
133    if (workvec3 == NULL)
134      {
135        printf("Storage Allocation Failed!\n");
136        exit(205);
137      };
138 
139    if (n > k)
140      {
141        workvec4=(double *)malloc(sizeof(double)*(n+1));
142      }
143    else
144      {
145        workvec4=(double *)malloc(sizeof(double)*(k+1));
146      };
147    if (workvec4 == NULL)
148      {
149        printf("Storage Allocation Failed!\n");
150        exit(205);
151      };
152 
153    if (n > k)
154      {
155        workvec5=(double *)malloc(sizeof(double)*(n+1));
156      }
157    else
158      {
159        workvec5=(double *)malloc(sizeof(double)*(k+1));
160      };
161    if (workvec5 == NULL)
162      {
163        printf("Storage Allocation Failed!\n");
164        exit(205);
165      };
166 
167    if (n > k)
168      {
169        workvec6=(double *)malloc(sizeof(double)*(n+1));
170      }
171    else
172      {
173        workvec6=(double *)malloc(sizeof(double)*(k+1));
174      };
175    if (workvec6 == NULL)
176      {
177        printf("Storage Allocation Failed!\n");
178        exit(205);
179      };
180 
181    if (n > k)
182      {
183        workvec7=(double *)malloc(sizeof(double)*(n+1));
184      }
185    else
186      {
187        workvec7=(double *)malloc(sizeof(double)*(k+1));
188      };
189    if (workvec7 == NULL)
190      {
191        printf("Storage Allocation Failed!\n");
192        exit(205);
193      };
194 
195    if (n > k)
196      {
197        workvec8=(double *)malloc(sizeof(double)*(n+1));
198      }
199    else
200      {
201        workvec8=(double *)malloc(sizeof(double)*(k+1));
202      };
203    if (workvec8 == NULL)
204      {
205        printf("Storage Allocation Failed!\n");
206        exit(205);
207      };
208 
209 
210    if (n > k)
211      {
212        diagO=(double *)malloc(sizeof(double)*(n+1));
213      }
214    else
215      {
216        diagO=(double *)malloc(sizeof(double)*(k+1));
217      };
218    if (diagO == NULL)
219      {
220        printf("Storage Allocation Failed!\n");
221        exit(205);
222      };
223 
224 
225 
226    rhs=malloc(sizeof(double)*(k+1));
227    if (rhs == NULL)
228      {
229        printf("Storage Allocation Failed!\n");
230        exit(205);
231      };
232 
233    dy=malloc(sizeof(double)*(k+1));
234    if (dy == NULL)
235      {
236        printf("Storage Allocation Failed!\n");
237        exit(205);
238      };
239 
240    dy1=malloc(sizeof(double)*(k+1));
241    if (dy1 == NULL)
242      {
243        printf("Storage Allocation Failed!\n");
244        exit(205);
245      };
246 
247    Fp=malloc(sizeof(double)*(k+1));
248    if (Fp == NULL)
249      {
250        printf("Storage Allocation Failed!\n");
251        exit(205);
252      };
253 
254    /*
255     * Check to make sure that there is at least one constraint.
256     */
257 
258    if (k<= 0)
259      {
260        printf("Problem must have at least one constraint.\n");
261        exit(206);
262      };
263 
264    /*
265     *  Work out the leading dimension for the array.  Note that we may not
266     *  want to use k itself, for cache issues.
267     */
268    if ((k % 2) == 0)
269      ldam=k+1;
270    else
271      ldam=k;
272 
273    O=malloc(sizeof(double)*ldam*ldam);
274    if (O == NULL)
275      {
276        printf("Storage Allocation Failed!\n");
277        exit(205);
278      };
279 
280    alloc_mat(C,&Zi);
281    alloc_mat(C,&dZ);
282    alloc_mat(C,&dX);
283 
284    /*
285     *  Fill in lots of details in the constraints data structure that haven't
286     *  necessarily been done before now.
287     */
288 
289    /*
290     * Next, setup issparse and NULL out all nextbyblock pointers.
291     */
292 
293    for (i=1; i<=k; i++)
294      {
295        p=constraints[i].blocks;
296        while (p != NULL)
297 	 {
298 	   /*
299 	    * First, set issparse.  The 0.125 parameter was hand tuned.
300 	    */
301 
302 	   if (((1.0*k*k*p->numentries*p->numentries) > 0.125*(1.0*k*p->blocksize*p->blocksize*p->blocksize)) && ((p->numentries) > 5))
303 	     {
304 	       p->issparse=0;
305 	     }
306 	   else
307 	     {
308 	       p->issparse=1;
309 	     };
310 
311 	   if (C.blocks[p->blocknum].blockcategory == DIAG)
312 	     p->issparse=1;
313 
314 	   /*
315 	    * Setup the cross links.
316 	    */
317 
318 	   p->nextbyblock=NULL;
319 	   p=p->next;
320 	 };
321      };
322 
323    /*
324     * Now, cross link.
325     */
326    for (i=1; i<=k; i++)
327      {
328        p=constraints[i].blocks;
329        while (p != NULL)
330 	 {
331 	   if (p->nextbyblock == NULL)
332 	     {
333 	       /*
334 		* link in the remaining blocks.
335 		*/
336 	       for (j=i+1; j<=k; j++)
337 		 {
338 		   q=constraints[j].blocks;
339 
340 		   while (q != NULL)
341 		     {
342 		       if (q->blocknum == p->blocknum)
343 			 {
344 			   if (p->nextbyblock == NULL)
345 			     {
346 			       p->nextbyblock=q;
347 			       q->nextbyblock=NULL;
348 			       prev=q;
349 			     }
350 			   else
351 			     {
352 			       prev->nextbyblock=q;
353 			       q->nextbyblock=NULL;
354 			       prev=q;
355 			     };
356 			   break;
357 			 };
358 		       q=q->next;
359 		     };
360 		 };
361 	     };
362 	   p=p->next;
363 	 };
364      };
365 
366    /*
367     * If necessary, print out information on sparsity of blocks.
368     */
369 
370    if (printlevel >= 1)
371      {
372        numblocks=0;
373        denseblocks=0;
374        for (i=1; i<=k; i++)
375 	 {
376 	   p=constraints[i].blocks;
377 	   while (p != NULL)
378 	     {
379                if (printlevel >= 4)
380                  printf("Sparsity of constraint blocks: %d,%d,%d,%d \n",i,p->blocknum,p->issparse,p->numentries);
381 
382                if (p->issparse == 0)
383                  denseblocks=denseblocks+1;
384 
385                numblocks=numblocks+1;
386 	       p=p->next;
387 	     };
388 	 };
389 
390        if (printlevel >= 2)
391          printf("Percentage of dense constraint blocks is %f\n",
392                 (100.0*denseblocks)/(1.0*numblocks));
393      };
394 
395    /*
396     * Allocate space for byblocks pointers.
397     */
398 
399    byblocks=(struct sparseblock **)malloc((C.nblocks+1)*sizeof(struct sparseblock *));
400    if (byblocks == NULL)
401      {
402        printf("Storage Allocation Failed!\n");
403        exit(205);
404      };
405 
406    for (i=1; i<=C.nblocks; i++)
407      byblocks[i]=NULL;
408 
409    /*
410     * Fill in byblocks pointers.
411     */
412 
413    for (i=1; i<=k; i++)
414      {
415        ptr=constraints[i].blocks;
416        while (ptr != NULL)
417 	 {
418 	   if (byblocks[ptr->blocknum]==NULL)
419 	     {
420 	       byblocks[ptr->blocknum]=ptr;
421 	     };
422 	   ptr=ptr->next;
423 	 };
424      };
425 
426    /*
427     *  Compute "fill".  This data structure tells us which elements in the
428     *  block diagonal matrix have corresponding elements in one of the
429     *  constraints, and which constraint this element first appears in.
430     *
431     */
432 
433    makefill(k,C,constraints,&fill,work1,printlevel);
434 
435    /*
436     * Compute the nonzero structure of O. This is only used for debugging
437     * at this point.
438     */
439 
440    if (printlevel >= 2)
441      {
442        nnz=structnnz(n,k,C,constraints);
443        printf("Structural density of O %d, %e \n",nnz,nnz*1.0/(k*k*1.0));
444      };
445 
446 
447    /*
448     * Sort entries in the constraints.
449     */
450 
451    sort_entries(k,C,constraints);
452 
453    /*
454     * Check the symmetry of C.
455     */
456 
457    checkc(n,C,printlevel);
458 
459    /*
460     * Check constraints.
461     */
462 
463    checkconstraints(n,k,C,constraints,printlevel);
464 
465    /*
466     *  Now, call sdp().
467     */
468 
469    ret=sdp(n,k,C,a,constant_offset,constraints,byblocks,fill,*pX,*py,*pZ,
470 	   cholxinv,cholzinv,ppobj,pdobj,work1,work2,work3,workvec1,
471 	   workvec2,workvec3,workvec4,workvec5,workvec6,workvec7,workvec8,
472 	   diagO,bestx,besty,bestz,Zi,O,rhs,dZ,dX,dy,dy1,Fp,
473 	   printlevel,params);
474 
475    if (printlevel >= 1)
476      {
477        if (ret==0)
478 	 printf("Success: SDP solved\n");
479        if (ret==1)
480 	 printf("Success: SDP is primal infeasible\n");
481        if (ret==2)
482 	 printf("Success: SDP is dual infeasible\n");
483        if (ret==3)
484 	 printf("Partial Success: SDP solved with reduced accuracy\n");
485        if (ret >= 4)
486 	 printf("Failure: return code is %d \n",ret);
487 
488        if (ret==1)
489 	 {
490 	   op_at(k,*py,constraints,work1);
491 	   addscaledmat(work1,-1.0,*pZ,work1);
492 	   printf("Certificate of primal infeasibility: a'*y=%.5e, ||A'(y)-Z||=%.5e\n",-1.0,Fnorm(work1));
493 	 };
494 
495        if (ret==2)
496 	 {
497 	   op_a(k,constraints,*pX,workvec1);
498 	   printf("Certificate of dual infeasibility: tr(CX)=%.5e, ||A(X)||=%.5e\n",trace_prod(C,*pX),norm2(k,workvec1+1));
499 	 };
500 
501        if ((ret==0) || (ret>=3))
502 	 {
503 	   if (printlevel >= 3)
504 	     {
505 	       printf("XZ Gap: %.7e \n",trace_prod(*pZ,*pX));
506 	       gap=*pdobj-*ppobj;
507 	       printf("Real Gap: %.7e \n",gap);
508 	     };
509 
510 	   if (printlevel >= 1)
511 	     {
512                gap=*pdobj-*ppobj;
513 
514 	       printf("Primal objective value: %.7e \n",*ppobj);
515 	       printf("Dual objective value: %.7e \n",*pdobj);
516 	       printf("Relative primal infeasibility: %.2e \n",
517 		      pinfeas(k,constraints,*pX,a,workvec1));
518 	       printf("Relative dual infeasibility: %.2e \n",
519 		      dinfeas(k,C,constraints,*py,*pZ,work1));
520 	       printf("Real Relative Gap: %.2e \n",gap/(1+fabs(*pdobj)+fabs(*ppobj)));
521 	       printf("XZ Relative Gap: %.2e \n",trace_prod(*pZ,*pX)/(1+fabs(*pdobj)+fabs(*ppobj)));
522 
523 	       printf("DIMACS error measures: %.2e %.2e %.2e %.2e %.2e %.2e\n",
524 		      pinfeas(k,constraints,*pX,a,workvec1)*(1+norm2(k,a+1))/
525 		      (1+norminf(k,a+1)),
526 		      0.0,
527 		      dimacserr3(k,C,constraints,*py,*pZ,work1),
528 		      0.0,
529 		      gap/(1+fabs(*pdobj)+fabs(*ppobj)),
530 		      trace_prod(*pZ,*pX)/(1+fabs(*pdobj)+fabs(*ppobj)));
531 	     };
532 	 };
533 
534      };
535 
536    /*
537     *  Now, free up all of the storage.
538     */
539 
540    free_mat(work1);
541    free_mat(work2);
542    free_mat(work3);
543    free_mat_packed(bestx);
544    free_mat_packed(bestz);
545    free_mat_packed(cholxinv);
546    free_mat_packed(cholzinv);
547 
548    free_mat(Zi);
549    free_mat(dZ);
550    free_mat(dX);
551 
552    free(besty);
553    free(workvec1);
554    free(workvec2);
555    free(workvec3);
556    free(workvec4);
557    free(workvec5);
558    free(workvec6);
559    free(workvec7);
560    free(workvec8);
561    free(rhs);
562    free(dy);
563    free(dy1);
564    free(Fp);
565    free(O);
566    free(diagO);
567    free(byblocks);
568 
569    /*
570     * Free up the fill data structure.
571     */
572 
573    ptr=fill.blocks;
574    while (ptr != NULL)
575      {
576        free(ptr->entries);
577        free(ptr->iindices);
578        free(ptr->jindices);
579        oldptr=ptr;
580        ptr=ptr->next;
581        free(oldptr);
582      };
583 
584 
585   /*
586    * Finally, free the constraints array.
587    */
588 
589    return(ret);
590 
591 }
592 
structnnz(n,k,C,constraints)593 int structnnz(n,k,C,constraints)
594      int n;
595      int k;
596      struct blockmatrix C;
597      struct constraintmatrix *constraints;
598 {
599   int i,j;
600   int ii,jj;
601   int nnz;
602   struct sparseblock *ptri;
603   struct sparseblock *ptrj;
604 
605 
606   nnz=0;
607   for (i=1; i<=k; i++)
608     for (j=1; j<=k; j++)
609       {
610 	ptri=constraints[i].blocks;
611 
612 	while (ptri != NULL)
613 	  {
614 
615 	    ptrj=constraints[j].blocks;
616 	    while (ptrj != NULL)
617 	      {
618 
619 		if (ptri->blocknum == ptrj->blocknum)
620 		  {
621 		    if (C.blocks[ptri->blocknum].blockcategory==MATRIX)
622 		      {
623 			nnz++;
624 			goto NEXTJ;
625 		      }
626 		    else
627 		      { /* DIAG block */
628 			for (ii=1; ii<=ptri->numentries; ii++)
629 			  for (jj=1; jj<=ptrj->numentries; jj++)
630 			    {
631 			      if (ptri->iindices[ii]==ptrj->iindices[jj])
632 				{
633 				  nnz++;
634 				  goto NEXTJ;
635 				};
636 			    };
637 		      };
638 		  };
639 		ptrj=ptrj->next;
640 	      };
641 
642 	    ptri=ptri->next;
643 	  }; /* end while */
644 
645       NEXTJ:;
646       }; /* end nested fors */
647 
648   return(nnz);
649 }
650 
actnnz(n,lda,A)651 int actnnz(n,lda,A)
652      int n;
653      int lda;
654      double A[];
655 {
656   int i,j;
657   int nnz;
658 
659   nnz=0;
660   for (i=1; i<=n; i++)
661     {
662       if (A[ijtok(i,i,lda)] != 0.0)
663 	nnz++;
664       for (j=i+1; j<=n; j++)
665 	{
666 	  if (A[ijtok(i,j,lda)] != 0.0)
667 	    {
668 	      nnz++;
669 	      nnz++;
670 	    };
671 	};
672     };
673 
674   return(nnz);
675 }
676 
bandwidth(n,lda,A)677 int bandwidth(n,lda,A)
678      int n;
679      int lda;
680      double A[];
681 {
682   int i;
683   int j;
684   int bw;
685 
686   bw=0;
687   for (j=2; j<=n; j++)
688     {
689       for (i=1; i<=j-1; i++)
690 	{
691 	  if (A[ijtok(i,j,lda)] != 0.0)
692 	    {
693 	      if ((j-i) > bw)
694 		bw=j-i;
695 	      break;
696 	    };
697 	};
698     };
699 
700   return(bw);
701 }
702 
703 /*
704  * Sanity checks for the C matrix data structure.  If this fails, we'll just
705  * exit(206) to indicate that there was improper input.  Otherwise, we'll
706  * return 0;
707  */
708 
checkc(int n,struct blockmatrix C,int printlevel)709 int checkc(int n,struct blockmatrix C,int printlevel)
710 {
711   int i,j,k;
712   int totalsize;
713 
714   totalsize=0;
715   for (k=1; k<=C.nblocks; k++)
716     {
717       if (C.blocks[k].blockcategory==DIAG)
718 	{
719 	  if (printlevel > 5)
720 	    printf("blockcategory=diag\n");
721 	};
722       if (C.blocks[k].blockcategory==MATRIX)
723 	{
724 	  if (printlevel > 5)
725 	    printf("blockcategory=matrix\n");
726 	};
727 
728       totalsize=totalsize+C.blocks[k].blocksize;
729 
730       if (C.blocks[k].blockcategory==MATRIX)
731 	{
732 	  for (i=1; i<=C.blocks[k].blocksize; i++)
733 	    {
734 	      for (j=1; j<=C.blocks[k].blocksize; j++)
735 		{
736 		  if (C.blocks[k].data.mat[ijtok(i,j,C.blocks[k].blocksize)] !=
737 		      C.blocks[k].data.mat[ijtok(j,i,C.blocks[k].blocksize)])
738 		    {
739                       if (printlevel >= 1)
740                         printf("C is not symmetric, %d, %d, %d\n",k,i,j);
741 		      exit(206);
742 		    }
743 		};
744 	    };
745 	};
746     };
747 
748   if (totalsize != n)
749     {
750       if (printlevel >= 1)
751         printf("Sum of block sizes does not equal n!\n");
752       exit(206);
753     };
754 
755   return(0);
756 }
757 
758 /*
759  * Sanity tests on the constraints data structure.
760  */
761 
checkconstraints(n,k,C,constraints,printlevel)762 int checkconstraints(n,k,C,constraints,printlevel)
763      int n;
764      int k;
765      struct blockmatrix C;
766      struct constraintmatrix *constraints;
767      int printlevel;
768 {
769   int i,j;
770   struct sparseblock *p;
771 
772   for (i=1; i<=k; i++)
773     {
774       p=constraints[i].blocks;
775       if (p==NULL)
776 	{
777           if (printlevel >= 1)
778             printf("Constraint %d is empty!\n",i);
779 	  exit(206);
780 	};
781       while (p != NULL)
782 	{
783 	  if (p->constraintnum != i)
784 	    {
785               if (printlevel >= 1)
786                 printf("p->constraintnum != i, i=%d \n",i);
787 	      exit(206);
788 	    };
789 	  if (p->blocksize != C.blocks[p->blocknum].blocksize)
790 	    {
791               if (printlevel >= 1)
792                 printf("p->blocksize is wrong, constraint %d \n",i);
793 	      exit(206);
794 	    };
795 	  if (printlevel > 5)
796 	    printf("Constraint %d, block %d, entries %d\n",i,p->blocknum,p->numentries);
797 	  for (j=1; j<=p->numentries; j++)
798 	    {
799 
800               /*
801                * For debugging, print out the constraint entry.
802                */
803 
804 	      if (printlevel >6)
805 		printf(" (%d, %d)=%lf\n",p->iindices[j],p->jindices[j],p->entries[j]);
806 
807               /*
808                * Make sure that all indices are in range.
809                */
810 
811 	      if (p->iindices[j] > C.blocks[p->blocknum].blocksize)
812 		{
813                   if (printlevel >= 1)
814                     printf("i index is larger than blocksize!\n");
815 		  exit(206);
816 		};
817 
818 	      if (p->jindices[j] > C.blocks[p->blocknum].blocksize)
819 		{
820                   if (printlevel >= 1)
821                     printf("j index is larger than blocksize!\n");
822 		  exit(206);
823 		};
824 
825 	      if (p->iindices[j] < 1)
826 		{
827                   if (printlevel >= 1)
828                     printf("i index is less than 1!\n");
829 		  exit(206);
830 		};
831 
832 	      if (p->jindices[j] < 1)
833 		{
834                   if (printlevel >= 1)
835                     printf("j index is less than 1!\n");
836 		  exit(206);
837 		};
838 
839               /*
840                * Make sure that entries are sorted properly.
841                */
842 
843               if (p->iindices[j] > p->jindices[j])
844                 {
845                   if (printlevel >= 1)
846                     {
847                       printf("i index is greater than j index!\n");
848                       printf("constraint=%d\n",i);
849                       printf("iindex=%d\n",p->iindices[j]);
850                       printf("jindex=%d\n",p->jindices[j]);
851                     };
852 
853 		  exit(206);
854                 };
855 
856 
857               /*
858                * Check for any duplicate entries that snuck in.
859                */
860 
861               if (j < p->numentries)
862                 if ((p->iindices[j]==p->iindices[j+1]) &
863                     (p->jindices[j]==p->jindices[j+1]))
864                   {
865                     if (printlevel >= 1)
866                       {
867                         printf("Duplicate entry!\n");
868                         printf("constraint=%d\n",i);
869                         printf("iindex=%d\n",p->iindices[j]);
870                         printf("jindex=%d\n",p->jindices[j]);
871                       };
872 
873                     exit(206);
874                   };
875 
876 	      if (C.blocks[p->blocknum].blockcategory==DIAG)
877 		{
878 		  if (p->iindices[j] != p->jindices[j])
879 		    {
880                       if (printlevel >= 1)
881                         printf("Off diagonal entry in diagonal block!\n");
882 		      exit(206);
883 		    };
884 		};
885 	    };
886 	  p=p->next;
887 	};
888     };
889   return(0);
890 }
891